1 Supplementary Material - Data checks, analysis and meta-a models

2 Data Set-up

2.1 Load packages

# Clear work space
  # rm(list=ls())

# Install CRAN packages
  library("pacman")
## Warning: package 'pacman' was built under R version 3.5.2
# Install orchard plot and metaAidR packages from GitHub
  devtools::install_github("itchyshin/orchard_plot", subdir = "orchaRd", force = TRUE, build_vignettes = TRUE); 
## Error in get(genname, envir = envir) : object 'testthat_print' not found
## 
##   
   checking for file ‘/private/var/folders/0b/pxghylq157gfhs1vrzdpx2gc0000gq/T/RtmpXG0rME/remotes6c851575cc49/itchyshin-orchard_plot-27d6281/orchaRd/DESCRIPTION’ ...
  
✓  checking for file ‘/private/var/folders/0b/pxghylq157gfhs1vrzdpx2gc0000gq/T/RtmpXG0rME/remotes6c851575cc49/itchyshin-orchard_plot-27d6281/orchaRd/DESCRIPTION’
## 
  
─  preparing ‘orchaRd’:
## 
  
   checking DESCRIPTION meta-information ...
  
✓  checking DESCRIPTION meta-information
## 
  
─  installing the package to build vignettes
## 
  
   creating vignettes ...
  
✓  creating vignettes (1.4s)
## 
  
─  checking for LF line-endings in source and make files and shell scripts
## 
  
─  checking for empty or unneeded directories
## ─  looking to see if a ‘data/datalist’ file should be added
## 
  
─  building ‘orchaRd_0.0.0.9000.tar.gz’
## 
  
   
## 
  devtools::install_github("daniel1noble/metaAidR"); 

pacman::p_load(knitr, metafor, dplyr, kableExtra, tidyverse, rotl, phytools, GGally, R.rsp, patchwork, devtools, robumeta, ape, geiger, phytools, phangorn, rlist, orchaRd, metaAidR, corrplot, stringr)

# set working directory
 # setwd("~/Documents/GitHub/sex_meta/")
  
# Source our own functions
  source("./R/func.R")
 
# Set the rerun object to FALSE so that you don't need to re-run all models again. Some take quite a lot of time to run. If FALSE, it will just re-load saved output. 

  rerun_models = FALSE

3 Organising data for analysis

  ## load our original pers dataset and our dataset with all of our sexual selection info
      pers <- read.csv("./data/pers_data.csv", stringsAsFactors = FALSE)
      bodysize <- read.csv("./data/bodysize_SSD.csv", stringsAsFactors = FALSE)

  ## Merge the two by spp_names columns
      pers <- merge(x = pers,
                   y = bodysize[,c("species_name", "SSD_index", "mating_system")],
                   by="species_name", all.x=TRUE,  no.dups = TRUE)

  ## Select only the relevant columns to make life easier
      pers_new <- pers %>% 
          select(study_ID, year, species_name, SSD_index, taxo_group, data_type, personality_trait, male_n, male_mean_conv, 
          male_sd_conv, female_n, female_mean_conv, female_sd_conv, depend, directionality, spp_name_phylo, mating_system, 
          age, population, study_environment, study_type, measurement_type)

  ## Add in observation level random effect (metafor doesn't do this, need to do it manually)
      pers_new <- pers_new %>% 
          group_by(taxo_group) %>% 
          mutate(obs = 1:length(study_ID))  

4 Calculating effect sizes

4.1 means - SMD using Hedges' g and variability - lnCVR

 ## SMD (Hedges' g)
   pers_new <- escalc(measure = "SMD", 
                     n1i = male_n, n2i = female_n,
                     m1i = male_mean_conv, m2i = female_mean_conv,
                     sd1i = male_sd_conv, sd2i = female_sd_conv, data = pers_new, var.names=c("SMD_yi","SMD_vi"), append = TRUE)
  
 ## lnCVR
   pers_new <- escalc(measure = "CVR",
                     n2i = female_n, n1i = male_n,
                     m2i = female_mean_conv, m1i = male_mean_conv,
                     sd2i = female_sd_conv, sd1i = male_sd_conv, data = pers_new, var.names=c("CVR_yi","CVR_vi"))

  # we have some NAs where one or both sexes have a value of 0 for either mean or SD. Will be easiest to just remove these.
    
   # Exclude NAs
      pers_new <- pers_new %>%
                    filter(!is.na(CVR_yi), !is.na(SMD_yi))
    
      dim(pers_new)    # check they've been removed with no issues

5 Data Checks

5.1 mean-variance relationship in our dataset

Looking at the strength of the correlation between the mean and SD to check that using lnCVR as a measure of variability is valid.
If the mean and SD are NOT strongly correlated then using lnCVR is pointless.

# females and males seperately because they are in different columns
   
   # use ggplot to make a scatterplot of females
    fem <- ggplot(pers_new, aes(x = female_mean_conv, y = female_sd_conv)) + geom_point()
    
    # on log scale
    fem + scale_x_continuous(trans = 'log10') + scale_y_continuous(trans = 'log10')

    # mean and SD on log scale to calculate correlation
    logfemale_mean <- log(pers_new$female_mean_conv)
    logfemale_SD <- log(pers_new$female_sd_conv)
    
    # correlation between mean and SD
    cor(logfemale_mean, logfemale_SD) #0.91
## [1] 0.9182668
  # Males
    # use ggplot to make a scatterplot of females
    male <- ggplot(pers_new, aes(x = male_mean_conv, y = male_sd_conv)) + geom_point()
    
    # on log scale
    male + scale_x_continuous(trans = 'log10') + scale_y_continuous(trans = 'log10')

    # mean and SD on log scale to calculate correlation
    logmale_mean <- log(pers_new$male_mean_conv)
    logmale_SD <- log(pers_new$male_sd_conv)
    
    # correlation between mean and SD
    cor(logmale_mean, logmale_SD) #0.90
## [1] 0.9071225

5.2 Checking for outliers and removing weird effect sizes

This is an important data checking step - here we can identify whether data has been entered or reported incorrectly (i.e. outliers)

First, let's look at the funnel plots for lnCVR and SMD. NOTE: these funnels have had our 2 big outliers already removed.

#funnel plot for lnCVR

    funnel(x = pers_new$CVR_yi, vi = pers_new$CVR_vi, yaxis="seinv", xlim = c(-10, 10))

SMD

#funnel plot for SMD

    funnel(x = pers_new$SMD_yi, vi = pers_new$SMD_vi, yaxis="seinv")
    text(as.character(pers_new$study_ID), x = pers_new$SMD_yi, y = 1/sqrt(pers_new$SMD_vi))

Removing outliers

    pers_new %>% 
    filter(study_ID == "P015" & SMD_yi < -15) # P015 has 1 large effect size, remove?
    
      # filter out that large effect size
      pers_new <- pers_new %>%
      filter(!study_ID == "P015" | !obs == "509")
      
      dim(pers_new) 
      
    # checking SMD outliers - inverse SE > 14
    funnel(x = pers_new$SMD_yi, vi = pers_new$SMD_vi, yaxis="seinv")
    text(as.character(pers_new$obs), x = pers_new$SMD_yi, y = 1/sqrt(pers_new$SMD_vi), offset = 0.8)
    
    # checking lnCVR outliers
    funnel(x = pers_new$CVR_yi, vi = pers_new$CVR_vi, yaxis="seinv", xlim = c(-10, 10))
    text(as.character(pers_new$obs), x = pers_new$CVR_yi, y = 1/sqrt(pers_new$CVR_vi), offset = 0.8) 
    
# Some measures are more physiological/not personality than personality, so probably wise to remove these before we run the models:
    # P029 - obs 22, 23, 32 
    # P084 - obs 59, 62, 63, 65, 68, 70, 71, 72, 74
    # P060 - obs 216, 217
    # P211 - obs 230, 245
    # P117 - obs 397, 393, 402, 414
    # P197 - obs 541, 544, 546, 547
    # P069 - obs 669, 672, 673, 682, 683, 684, 686, 694 
    
    # remove these by study 
    pers_new <- pers_new %>% filter(!study_ID == "P029" | !obs %in% c("21", "25", "26", "28", "31"))
    
    pers_new <- pers_new %>% filter(!study_ID == "P084" | !obs %in% c("59", "62", "63", "65", "68", "70", "71", "72", "74"))
    
    pers_new <-  pers_new %>% filter(!study_ID == "P060" | !obs %in% c("216", "217"))
    
    pers_new <- pers_new %>% filter(!study_ID == "P211" | !obs %in% c("230", "245"))
    
    pers_new <- pers_new %>% filter(!study_ID == "P117" | !obs %in% c("397", "393", "402", "414"))
    
    pers_new <- pers_new %>% filter(!study_ID == "P197" | !obs %in% c("541", "544", "546", "547"))
    
    pers_new <- pers_new %>% filter(!study_ID == "P197" | !obs %in% c("669", "672", "673", "682", "683", "684", "686", "694")) 
    pers_new <- pers_new %>% filter(!study_ID == "P041" | !obs %in% c("120", "124"))
    
    dim(pers_new) # check they've been removed without issue
    
    # after checking the data, there are a few effect sizes that might be driving weird results so let's drop them and see
    
    pers_new <- pers_new %>% filter(!study_ID == "P100" | !obs == "519") # big outlier

5.3 Flip signs of effects for SMD

The directional meaning of effect sizes vary depending on the specific units and trait being measured. The data has a directionality column that tells one if the meaning should be reversed (1) or left the same.

   pers_new$directionality <- ifelse(is.na(pers_new$directionality), 0, 1)

      pers_new$SMD_yi_flip <- ifelse(pers_new$directionality == 1, pers_new$SMD_yi*(-1), pers_new$SMD_yi)

6 Prepare the phylogenetic trees

We constructed seperate phylogenetic trees for each taxonomic group. The tree for birds was constructed using BirdTree.org, the rest were constructed using TimeTree.org. We'll use these trees for multi-level meta-analytic models throughout the analysis.

# Find all tree file names
  tree_files <- paste0("./trees/", list.files("./trees"))[-1]
  
    # Bird tree has been constructed already, just need to get trees for the rest of the taxo groups   
      trees <- lapply(tree_files, function(x) read.tree(x))
      names <- gsub("~./trees/", "", tree_files)
      names(trees) <- names

    # Plot the trees and see how they look
      par(mfrow = c(1,5), mar = c(1,1,1,1))
      lapply(trees, function(x) plot(x, cex = 1))

## $`./trees/bird_species.nwk`
## $`./trees/bird_species.nwk`$type
## [1] "phylogram"
## 
## $`./trees/bird_species.nwk`$use.edge.length
## [1] TRUE
## 
## $`./trees/bird_species.nwk`$node.pos
## [1] 1
## 
## $`./trees/bird_species.nwk`$node.depth
## [1] 1
## 
## $`./trees/bird_species.nwk`$show.tip.label
## [1] TRUE
## 
## $`./trees/bird_species.nwk`$show.node.label
## [1] FALSE
## 
## $`./trees/bird_species.nwk`$font
## [1] 3
## 
## $`./trees/bird_species.nwk`$cex
## [1] 1
## 
## $`./trees/bird_species.nwk`$adj
## [1] 0
## 
## $`./trees/bird_species.nwk`$srt
## [1] 0
## 
## $`./trees/bird_species.nwk`$no.margin
## [1] FALSE
## 
## $`./trees/bird_species.nwk`$label.offset
## [1] 0
## 
## $`./trees/bird_species.nwk`$x.lim
## [1]    0.000 1716.674
## 
## $`./trees/bird_species.nwk`$y.lim
## [1]   1 109
## 
## $`./trees/bird_species.nwk`$direction
## [1] "rightwards"
## 
## $`./trees/bird_species.nwk`$tip.color
## [1] "black"
## 
## $`./trees/bird_species.nwk`$Ntip
## [1] 109
## 
## $`./trees/bird_species.nwk`$Nnode
## [1] 108
## 
## $`./trees/bird_species.nwk`$root.time
## NULL
## 
## $`./trees/bird_species.nwk`$align.tip.label
## [1] FALSE
## 
## 
## $`./trees/fish_species.nwk`
## $`./trees/fish_species.nwk`$type
## [1] "phylogram"
## 
## $`./trees/fish_species.nwk`$use.edge.length
## [1] TRUE
## 
## $`./trees/fish_species.nwk`$node.pos
## [1] 1
## 
## $`./trees/fish_species.nwk`$node.depth
## [1] 1
## 
## $`./trees/fish_species.nwk`$show.tip.label
## [1] TRUE
## 
## $`./trees/fish_species.nwk`$show.node.label
## [1] FALSE
## 
## $`./trees/fish_species.nwk`$font
## [1] 3
## 
## $`./trees/fish_species.nwk`$cex
## [1] 1
## 
## $`./trees/fish_species.nwk`$adj
## [1] 0
## 
## $`./trees/fish_species.nwk`$srt
## [1] 0
## 
## $`./trees/fish_species.nwk`$no.margin
## [1] FALSE
## 
## $`./trees/fish_species.nwk`$label.offset
## [1] 0
## 
## $`./trees/fish_species.nwk`$x.lim
## [1]     0.00 35169.86
## 
## $`./trees/fish_species.nwk`$y.lim
## [1]  1 22
## 
## $`./trees/fish_species.nwk`$direction
## [1] "rightwards"
## 
## $`./trees/fish_species.nwk`$tip.color
## [1] "black"
## 
## $`./trees/fish_species.nwk`$Ntip
## [1] 22
## 
## $`./trees/fish_species.nwk`$Nnode
## [1] 21
## 
## $`./trees/fish_species.nwk`$root.time
## NULL
## 
## $`./trees/fish_species.nwk`$align.tip.label
## [1] FALSE
## 
## 
## $`./trees/invert_species.nwk`
## $`./trees/invert_species.nwk`$type
## [1] "phylogram"
## 
## $`./trees/invert_species.nwk`$use.edge.length
## [1] TRUE
## 
## $`./trees/invert_species.nwk`$node.pos
## [1] 1
## 
## $`./trees/invert_species.nwk`$node.depth
## [1] 1
## 
## $`./trees/invert_species.nwk`$show.tip.label
## [1] TRUE
## 
## $`./trees/invert_species.nwk`$show.node.label
## [1] FALSE
## 
## $`./trees/invert_species.nwk`$font
## [1] 3
## 
## $`./trees/invert_species.nwk`$cex
## [1] 1
## 
## $`./trees/invert_species.nwk`$adj
## [1] 0
## 
## $`./trees/invert_species.nwk`$srt
## [1] 0
## 
## $`./trees/invert_species.nwk`$no.margin
## [1] FALSE
## 
## $`./trees/invert_species.nwk`$label.offset
## [1] 0
## 
## $`./trees/invert_species.nwk`$x.lim
## [1]    0.000 8492.817
## 
## $`./trees/invert_species.nwk`$y.lim
## [1]  1 42
## 
## $`./trees/invert_species.nwk`$direction
## [1] "rightwards"
## 
## $`./trees/invert_species.nwk`$tip.color
## [1] "black"
## 
## $`./trees/invert_species.nwk`$Ntip
## [1] 42
## 
## $`./trees/invert_species.nwk`$Nnode
## [1] 41
## 
## $`./trees/invert_species.nwk`$root.time
## NULL
## 
## $`./trees/invert_species.nwk`$align.tip.label
## [1] FALSE
## 
## 
## $`./trees/mammal_species.nwk`
## $`./trees/mammal_species.nwk`$type
## [1] "phylogram"
## 
## $`./trees/mammal_species.nwk`$use.edge.length
## [1] TRUE
## 
## $`./trees/mammal_species.nwk`$node.pos
## [1] 1
## 
## $`./trees/mammal_species.nwk`$node.depth
## [1] 1
## 
## $`./trees/mammal_species.nwk`$show.tip.label
## [1] TRUE
## 
## $`./trees/mammal_species.nwk`$show.node.label
## [1] FALSE
## 
## $`./trees/mammal_species.nwk`$font
## [1] 3
## 
## $`./trees/mammal_species.nwk`$cex
## [1] 1
## 
## $`./trees/mammal_species.nwk`$adj
## [1] 0
## 
## $`./trees/mammal_species.nwk`$srt
## [1] 0
## 
## $`./trees/mammal_species.nwk`$no.margin
## [1] FALSE
## 
## $`./trees/mammal_species.nwk`$label.offset
## [1] 0
## 
## $`./trees/mammal_species.nwk`$x.lim
## [1]    0.000 1711.683
## 
## $`./trees/mammal_species.nwk`$y.lim
## [1]  1 46
## 
## $`./trees/mammal_species.nwk`$direction
## [1] "rightwards"
## 
## $`./trees/mammal_species.nwk`$tip.color
## [1] "black"
## 
## $`./trees/mammal_species.nwk`$Ntip
## [1] 46
## 
## $`./trees/mammal_species.nwk`$Nnode
## [1] 45
## 
## $`./trees/mammal_species.nwk`$root.time
## NULL
## 
## $`./trees/mammal_species.nwk`$align.tip.label
## [1] FALSE
## 
## 
## $`./trees/reptile_species.nwk`
## $`./trees/reptile_species.nwk`$type
## [1] "phylogram"
## 
## $`./trees/reptile_species.nwk`$use.edge.length
## [1] TRUE
## 
## $`./trees/reptile_species.nwk`$node.pos
## [1] 1
## 
## $`./trees/reptile_species.nwk`$node.depth
## [1] 1
## 
## $`./trees/reptile_species.nwk`$show.tip.label
## [1] TRUE
## 
## $`./trees/reptile_species.nwk`$show.node.label
## [1] FALSE
## 
## $`./trees/reptile_species.nwk`$font
## [1] 3
## 
## $`./trees/reptile_species.nwk`$cex
## [1] 1
## 
## $`./trees/reptile_species.nwk`$adj
## [1] 0
## 
## $`./trees/reptile_species.nwk`$srt
## [1] 0
## 
## $`./trees/reptile_species.nwk`$no.margin
## [1] FALSE
## 
## $`./trees/reptile_species.nwk`$label.offset
## [1] 0
## 
## $`./trees/reptile_species.nwk`$x.lim
## [1]    0.000 2656.912
## 
## $`./trees/reptile_species.nwk`$y.lim
## [1]  1 10
## 
## $`./trees/reptile_species.nwk`$direction
## [1] "rightwards"
## 
## $`./trees/reptile_species.nwk`$tip.color
## [1] "black"
## 
## $`./trees/reptile_species.nwk`$Ntip
## [1] 10
## 
## $`./trees/reptile_species.nwk`$Nnode
## [1] 9
## 
## $`./trees/reptile_species.nwk`$root.time
## NULL
## 
## $`./trees/reptile_species.nwk`$align.tip.label
## [1] FALSE
# Checking trees to ensure we only include species in the current dataset

    # Check that they are ultrametric
      lapply(trees, function(x) is.ultrametric(x))

    # Check that all names in the phylogeny are also in the data
      taxa_data_list <- split(pers_new, pers_new$taxo_group)
      
      other_groups <- mapply(x = taxa_data_list, 
                             y = trees, 
                             function(x,y) tree_checks(x,y, "spp_name_phylo", type = "checks")) 

    # Print out each taxon group
      for(i in colnames(other_groups)){
            print(i)
            print(other_groups[,i] )
      }
      

      # Now to prune trees so that we get tree names that match with species in data
      pruned_trees <- mapply(x = taxa_data_list, 
                             y = trees, 
                             function(x,y) tree_checks(x,y, "spp_name_phylo", type = "prune"))

      # Check that this has been done correctly
        re_checks <- mapply(x = taxa_data_list, 
                             y = pruned_trees, 
                             function(x,y) tree_checks(x,y, "spp_name_phylo", type = "checks"))
        
      for(i in colnames(re_checks)){
            print(i)
            print(re_checks[,i] )
      }

# Extract the phylogenetic correlation matrices
      phylo_vcv <- lapply(pruned_trees, function(x) vcv(x, corr = TRUE)) # these matrices are used in the meta-a models

7 Sensitivity analysis - checking score data

Before we begin, we need to run a sensitivity analysis to see if score data is ok to use. With these models, we are just including score as a moderator term to compare with the rest of the dataset (some of which has already been transformed, we just can't do that with scores). Model summaries are also presented in Supplementary Table S2.

Our score sensitivity model:

    # model:
     sensitivity_mod1_score <- meta_model_fits(pers_new, phylo_vcv, type = "score")

    # Extract the SMD and lnCVR results
     smd_mods_score <- sensitivity_mod1_score["SMD",]  
     lnCVR_mods_score <- sensitivity_mod1_score["lnCVR",] # inverts significant
    
    # Because invert score data is significantly different, we need to remove these effect sizes before running our models
        
      # filter out invert scores from pers dataset  
            pers_new <- pers_new %>%
                              filter(score != "score" | taxo_group != "invertebrate") 
    
            dim(pers_new) 
            
   # How many effect sizes, unique studies and different species are we left with?
            data.frame(pers_new %>% 
                         summarise(n = n(), studies = length(unique(study_ID)), species = length(unique(spp_name_phylo))))

8 Meta-Analysis Models

Let's run the first bunch of models on the whole dataset. We'll start off with intercept-only multi-level meta-analytic models, then move to multi-level meta-regression models (personality traits, and SSD). The functions in func.R should be consulted to see precisely what models are being fit across the taxonomic groups.

8.1 Intercept-only MLMA models

Complete model summaries are also presented in Supplementary Table S14.

# First we will fit our MLMA intercept only models, across each taxo group. 

  # we can use this function to just read the saved model output instead of re-running the model, which takes a while
rerun_models == FALSE
## [1] TRUE
  if(rerun_models == TRUE){
    MLMA_models <- meta_model_fits(pers_new, phylo_vcv, type = "int")
    saveRDS(MLMA_models, "./output/MLMA_models_int")
  }else{
    MLMA_models <- readRDS("./output/MLMA_models_int")
  }


  # View model results
  split_taxa <- split(pers_new, pers_new$taxo_group)
  
    smd_mods <- MLMA_models["SMD",]
  
    lnCVR_mods <- MLMA_models["lnCVR",]

8.2 I2 estimates of heterogeneity - intercept models

Study_ID is the between study heterogeneity, Phylo tells us if there is a phylogenetic signal and the strength of that signal. Total I2 is testing how much heterogeneity we have beyond sampling variance

# From these models we can get I2 estimates:
  birds_smd = I2(smd_mods[[1]], v = split_taxa[[1]]$SMD_vi, phylo = "spp_name_phylo", obs = "obs")
     
  birds_CVR = I2(lnCVR_mods[[1]], v = split_taxa[[1]]$CVR_vi, phylo = "spp_name_phylo", obs = "obs")
      
  fish_smd = I2(smd_mods[[2]], v = split_taxa[[2]]$SMD_vi, phylo = "spp_name_phylo", obs = "obs")
      
  fish_CVR = I2(lnCVR_mods[[2]], v = split_taxa[[2]]$CVR_vi, phylo = "spp_name_phylo", obs = "obs")
      
  invert_smd = I2(smd_mods[[3]], v = split_taxa[[3]]$SMD_vi, phylo = "spp_name_phylo", obs = "obs")
      
  invert_CVR = I2(lnCVR_mods[[3]], v = split_taxa[[3]]$CVR_vi, phylo = "spp_name_phylo", obs = "obs")
      
  mammal_smd = I2(smd_mods[[4]], v = split_taxa[[4]]$SMD_vi, phylo = "spp_name_phylo", obs = "obs")
      
  mammal_CVR = I2(lnCVR_mods[[4]], v = split_taxa[[4]]$CVR_vi, phylo = "spp_name_phylo", obs = "obs")
      
  reptile_smd = I2(smd_mods[[5]], v = split_taxa[[5]]$SMD_vi, phylo = "spp_name_phylo", obs = "obs")
      
  reptile_CVR = I2(lnCVR_mods[[5]], v = split_taxa[[5]]$CVR_vi, phylo = "spp_name_phylo", obs = "obs")
      
  
# Now that we have our list of models, we can extract the estimates, CIs and prediction intervals
  MLMA_estimates_SMD <- plyr::ldply(lapply(smd_mods, function(x) print(mod_results(x, mod = "Int"))))
      MLMA_estimates_SMD
  MLMA_estimates_lnCVR <- plyr::ldply(lapply(lnCVR_mods, function(x) print(mod_results(x, mod = "Int"))))
      MLMA_estimates_lnCVR

Extract p-values from these models to use later when adjusting them for multiple testing:

# taking p-values from models for False Discovery Rate p-value adjustment
  p.SMD_intercept <- unlist(lapply(smd_mods, function(x) x$pval))
      p.SMD_intercept
##         bird         fish invertebrate       mammal     reptilia 
##    0.3560563    0.3900624    0.1778479    0.7143065    0.4793203
  p.lnCVR_intercept <- unlist(lapply(lnCVR_mods, function(x) x$pval))
      p.lnCVR_intercept
##         bird         fish invertebrate       mammal     reptilia 
##    0.5812731    0.9117446    0.8659517    0.6776659    0.3554154

8.3 Personality trait MLMR models

These models include personality trait type as a moderator. Please note that we estimate the mean for each of the categorical levels because we are not really interested in whether the means differ, but whether or not males and females differ in any of these traits.

Complete model summaries are presented in Supplementary Table S15.

# we can just reload saved model outputs here to save time
rerun_models == FALSE
## [1] TRUE
    if(rerun_models == TRUE){
      MLMR_models_pers_trait <- meta_model_fits(pers_new, phylo_vcv, type = "pers")
      saveRDS(MLMR_models_pers_trait, "./output/MLMR_models_pers_trait")
    } else{
      MLMR_models_pers_trait <- readRDS("./output/MLMR_models_pers_trait")
    }

  # Extract the SMD and lnCVR results
  smd_mods_pers <- MLMR_models_pers_trait["SMD",]
     
  lnCVR_mods_pers <- MLMR_models_pers_trait["lnCVR",] 
     
  # these model objects are used to make the orchard plots shown in Figures 2-6

Get prediction intervals for personality trait models:

  # Get the combined estimates from them all
  MLMA_estimates_SMD_pers <- plyr::ldply(lapply(smd_mods_pers, function(x) print(mod_results(x, mod = "personality_trait"))))
      MLMA_estimates_SMD_pers
  MLMA_estimates_lnCVR_pers <- plyr::ldply(lapply(lnCVR_mods_pers, function(x) print(mod_results(x, mod = "personality_trait"))))
      MLMA_estimates_lnCVR_pers
  # Add in n and k to these dataframes
  n_k<- pers_new %>%
      group_by(taxo_group, personality_trait) %>%
      summarise(n = n(), spp = length(unique(spp_name_phylo)), k = length(unique(study_ID)))
  
  # Summary of model estimates with number of studies, species and effect sizes included
  MLMA_estimates_SMD_pers <- data.frame(MLMA_estimates_SMD_pers, n_k[,c("n", "spp", "k")])
      MLMA_estimates_SMD_pers
  MLMA_estimates_lnCVR_pers <- data.frame(MLMA_estimates_lnCVR_pers, n_k[,c("n", "spp", "k")])
      MLMA_estimates_lnCVR_pers

Extract p-values from models for multiple testing adjustment later:

# extract p-values for multiple testing 
  p.SMD_pers <- unlist(lapply(smd_mods_pers, function(x) x$pval))
    p.SMD_pers
##         bird1         bird2         bird3         bird4         bird5 
##    0.57653039    0.30479363    0.27751708    0.65883812    0.01348632 
##         fish1         fish2         fish3         fish4         fish5 
##    0.51728620    0.27173344    0.70897900    0.26038453    0.54643814 
## invertebrate1 invertebrate2 invertebrate3 invertebrate4 invertebrate5 
##    0.21528997    0.23621794    0.22305017    0.86210097    0.20535070 
##       mammal1       mammal2       mammal3       mammal4       mammal5 
##    0.42949385    0.64539838    0.47941582    0.89056692    0.74385107 
##     reptilia1     reptilia2     reptilia3     reptilia4     reptilia5 
##    0.83127102    0.45179068    0.38700906    0.03843513    0.98694696
  p.lnCVR_pers <- unlist(lapply(lnCVR_mods_pers, function(x) x$pval))
    p.lnCVR_pers
##         bird1         bird2         bird3         bird4         bird5 
##    0.66827246    0.63819433    0.71794905    0.02062033    0.56829592 
##         fish1         fish2         fish3         fish4         fish5 
##    0.89295129    0.08859994    0.67087141    0.63106746    0.02038126 
## invertebrate1 invertebrate2 invertebrate3 invertebrate4 invertebrate5 
##    0.54416559    0.80512014    0.72746311    0.21252880    0.56927261 
##       mammal1       mammal2       mammal3       mammal4       mammal5 
##    0.46265634    0.45330729    0.82211806    0.72515472    0.76995584 
##     reptilia1     reptilia2     reptilia3     reptilia4     reptilia5 
##    0.49192601    0.06178802    0.20080617    0.25517952    0.90971983

8.4 Personality trait x SSD MLMR models

Now let's look at how SSD interacts with personality trait type. Here we are not estimating an intercept either, so each intercept varies by trait category and each slope as well. Note that there are lots of warnings, but these are the result of many levels not being present in taxa groups. We chose not to scale SSD_index because it is easier (and biologically relevant) to interpret SSD when it is 0 (when males and females are the same size), and when SSD is positive (when males are larger than females).

Model summaries are presented in Supplementary Table S17.

# again, we can just reload our saved model output here  
rerun_models == FALSE
## [1] TRUE
      if(rerun_models == TRUE){
      MLMR_models_pers_SSD <- meta_model_fits(pers_new, phylo_vcv, type = "pers_SSD")
      saveRDS(MLMR_models_pers_SSD, "./output/MLMR_models_pers_SSD")
    } else{
      MLMR_models_pers_SSD <- readRDS("./output/MLMR_models_pers_SSD")
    }  


  # Extract the SMD and lnCVR results
  smd_mods_pers_SSD <- MLMR_models_pers_SSD["SMD",]
    
  lnCVR_mods_pers_SSD <- MLMR_models_pers_SSD["lnCVR",]

Get the prediction intervals for our interaction models:

  # extract estimates using modified function in func.R file:
    # SMD
  MLMA_estimates_SMD_SSD <- plyr::ldply(lapply(smd_mods_pers_SSD, function(x) 
     print(mod_results_new(x, mod_cat = "personality_trait", mod_cont = "SSD_index", type = "zero"))))
    
    MLMA_estimates_SMD_SSD
    # lnCVR
  MLMA_estimates_lnCVR_SSD <- plyr::ldply(lapply(lnCVR_mods_pers_SSD, function(x) 
   print(mod_results_new(x, mod_cat = "personality_trait", mod_cont = "SSD_index", type = "zero"))))
    
    MLMA_estimates_lnCVR_SSD
  # Table to get species numbers, no. studies and no. effect sizes:
  data.frame(pers_new %>%
  group_by(taxo_group, personality_trait) %>%
  filter(!is.na(SSD_index))%>%
  summarise(n = n(), N_spp = length(unique(spp_name_phylo)), N_studies = length(unique(study_ID))))

8.5 SSD subset models

Because we aren't really interested in how each trait type differs from each other, we need to run our SSD models on subsets of the data where we can get the mean estimates for individual trait types and for SSD. Model summaries are presented in Supplementary Table S16.

NOTE: Since we are conducting our meta-regression at the species level (the level at which we can assume effect sizes are independent), any personality trait with fewer than 10 species needs to be dropped to look at interactions between SSD and personality. Having a minimum of 10 studies etc. is the rule of thumb for meta-regressions (e.g. see Borenstein et al Intro to Meta-A)

8.5.1 NEED TO DROP:

  1. ALL REPTILES - GONE, NOT ENOUGH SPECIES
  2. BIRDS - EVERYTHING BUT BOLDNESS
  3. FISH - ACTIVITY, EXPLORATION & SOCIALITY
  4. INVERTS - SOCIALITY, EXPLORATION & AGGRESSION
  5. MAMMALS - SOCIALITY

8.5.2 Mammals

SSD for activity, boldness, aggression and exploration.

# 1. MAMMALS

# First, we need to subset our pers dataset by taxo group to drop the unwanted levels.
  # a. activity
    pers_new_mammal_activity <- as.data.frame(pers_new %>%
    filter(personality_trait == "activity") %>%
    filter(taxo_group == "mammal")) 
  # b. boldness
    pers_new_mammal_boldness <- as.data.frame(pers_new %>%
    filter(personality_trait == "boldness") %>%
    filter(taxo_group == "mammal")) 
  # c. aggression
    pers_new_mammal_aggression <- as.data.frame(pers_new %>%
    filter(personality_trait == "aggression") %>%
    filter(taxo_group == "mammal")) 
  # d. exploration
    pers_new_mammal_exploration <- as.data.frame(pers_new %>%
    filter(personality_trait == "exploration") %>%
    filter(taxo_group == "mammal")) 

  # Extract the phylogenetic correlation matrices
     phylo_vcv_mammal <- phylo_vcv[[4]] 

Activity:

    # a. activity
    #SMD
    MLMR_mods_pers_SSD_mammal_activity_SMD <- rma.mv(SMD_yi_flip ~ SSD_index, V = SMD_vi, 
                                          random = list(~1|study_ID, ~1|spp_name_phylo, ~1|obs), 
                                          R = list(spp_name_phylo=phylo_vcv_mammal), control=list(optimizer="optim"), 
                                          test = "t", data = pers_new_mammal_activity)
## Warning in rma.mv(SMD_yi_flip ~ SSD_index, V = SMD_vi, random = list(~1 | :
## There are rows/columns in the 'R' matrix for 'spp_name_phylo' for which there
## are no data.
    MLMR_mods_pers_SSD_mammal_activity_SMD
## 
## Multivariate Meta-Analysis Model (k = 83; method: REML)
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed          factor    R 
## sigma^2.1  0.1000  0.3163     14     no        study_ID   no 
## sigma^2.2  2.8075  1.6755     12     no  spp_name_phylo  yes 
## sigma^2.3  0.1914  0.4375     83     no             obs   no 
## 
## Test for Residual Heterogeneity:
## QE(df = 81) = 320.2151, p-val < .0001
## 
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 81) = 5.0697, p-val = 0.0271
## 
## Model Results:
## 
##            estimate      se     tval  df    pval    ci.lb    ci.ub 
## intrcpt      0.5040  1.2531   0.4022  81  0.6886  -1.9892   2.9972    
## SSD_index   -2.2054  0.9795  -2.2516  81  0.0271  -4.1543  -0.2565  * 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    #lnCVR
    MLMR_mods_pers_SSD_mammal_activity_lncvr <- rma.mv(CVR_yi ~ SSD_index, V = CVR_vi, 
                                            random = list(~1|study_ID, ~1|spp_name_phylo, ~1|obs), 
                                            R = list(spp_name_phylo=phylo_vcv_mammal), control=list(optimizer="optim"), 
                                            test = "t", data = pers_new_mammal_activity)
## Warning in rma.mv(CVR_yi ~ SSD_index, V = CVR_vi, random = list(~1 | study_ID, :
## There are rows/columns in the 'R' matrix for 'spp_name_phylo' for which there
## are no data.
    MLMR_mods_pers_SSD_mammal_activity_lncvr
## 
## Multivariate Meta-Analysis Model (k = 83; method: REML)
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed          factor    R 
## sigma^2.1  0.0293  0.1710     14     no        study_ID   no 
## sigma^2.2  0.0001  0.0079     12     no  spp_name_phylo  yes 
## sigma^2.3  0.0607  0.2465     83     no             obs   no 
## 
## Test for Residual Heterogeneity:
## QE(df = 81) = 146.2596, p-val < .0001
## 
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 81) = 0.1324, p-val = 0.7169
## 
## Model Results:
## 
##            estimate      se    tval  df    pval    ci.lb   ci.ub 
## intrcpt      0.0518  0.0988  0.5236  81  0.6020  -0.1449  0.2484    
## SSD_index    0.1248  0.3430  0.3638  81  0.7169  -0.5577  0.8073    
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Boldness:

    # b. boldness
    #SMD
    MLMR_mods_pers_SSD_mammal_bold_SMD <- rma.mv(SMD_yi_flip ~ SSD_index, V = SMD_vi, 
                                              random = list(~1|study_ID, ~1|spp_name_phylo, ~1|obs), 
                                              R = list(spp_name_phylo=phylo_vcv_mammal), control=list(optimizer="optim"), 
                                              test = "t", data = pers_new_mammal_boldness)
## Warning in rma.mv(SMD_yi_flip ~ SSD_index, V = SMD_vi, random = list(~1 | :
## There are rows/columns in the 'R' matrix for 'spp_name_phylo' for which there
## are no data.
## Warning: Rows with NAs omitted from model fitting.
    MLMR_mods_pers_SSD_mammal_bold_SMD    
## 
## Multivariate Meta-Analysis Model (k = 163; method: REML)
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed          factor    R 
## sigma^2.1  0.0088  0.0938     26     no        study_ID   no 
## sigma^2.2  0.0000  0.0050     26     no  spp_name_phylo  yes 
## sigma^2.3  0.1707  0.4132    163     no             obs   no 
## 
## Test for Residual Heterogeneity:
## QE(df = 161) = 405.7659, p-val < .0001
## 
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 161) = 1.3101, p-val = 0.2541
## 
## Model Results:
## 
##            estimate      se     tval   df    pval    ci.lb   ci.ub 
## intrcpt      0.0739  0.0774   0.9547  161  0.3412  -0.0789  0.2267    
## SSD_index   -0.1686  0.1473  -1.1446  161  0.2541  -0.4594  0.1223    
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    #lnCVR
    MLMR_mods_pers_SSD_mammal_bold_lncvr <- rma.mv(CVR_yi ~ SSD_index, V = CVR_vi, 
                                            random = list(~1|study_ID, ~1|spp_name_phylo, ~1|obs), 
                                            R = list(spp_name_phylo=phylo_vcv_mammal), control=list(optimizer="optim"), 
                                            test = "t", data = pers_new_mammal_boldness)
## Warning in rma.mv(CVR_yi ~ SSD_index, V = CVR_vi, random = list(~1 | study_ID, :
## There are rows/columns in the 'R' matrix for 'spp_name_phylo' for which there
## are no data.

## Warning in rma.mv(CVR_yi ~ SSD_index, V = CVR_vi, random = list(~1 | study_ID, :
## Rows with NAs omitted from model fitting.
    MLMR_mods_pers_SSD_mammal_bold_lncvr
## 
## Multivariate Meta-Analysis Model (k = 163; method: REML)
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed          factor    R 
## sigma^2.1  0.0029  0.0543     26     no        study_ID   no 
## sigma^2.2  0.0000  0.0027     26     no  spp_name_phylo  yes 
## sigma^2.3  0.0211  0.1452    163     no             obs   no 
## 
## Test for Residual Heterogeneity:
## QE(df = 161) = 177.7499, p-val = 0.1737
## 
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 161) = 2.1066, p-val = 0.1486
## 
## Model Results:
## 
##            estimate      se    tval   df    pval    ci.lb   ci.ub 
## intrcpt      0.0114  0.0513  0.2230  161  0.8238  -0.0899  0.1128    
## SSD_index    0.1316  0.0907  1.4514  161  0.1486  -0.0475  0.3106    
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Aggression:

    # c. aggression
    #SMD
    MLMR_mods_pers_SSD_mammal_aggression_SMD <- rma.mv(SMD_yi_flip ~ SSD_index, V = SMD_vi, 
                                          random = list(~1|study_ID, ~1|spp_name_phylo, ~1|obs), 
                                          R = list(spp_name_phylo=phylo_vcv_mammal), control=list(optimizer="optim"), 
                                          test = "t", data = pers_new_mammal_aggression)
## Warning in rma.mv(SMD_yi_flip ~ SSD_index, V = SMD_vi, random = list(~1 | :
## There are rows/columns in the 'R' matrix for 'spp_name_phylo' for which there
## are no data.
## Warning: Rows with NAs omitted from model fitting.
    MLMR_mods_pers_SSD_mammal_aggression_SMD
## 
## Multivariate Meta-Analysis Model (k = 85; method: REML)
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed          factor    R 
## sigma^2.1  0.0000  0.0039     15     no        study_ID   no 
## sigma^2.2  0.6852  0.8277     13     no  spp_name_phylo  yes 
## sigma^2.3  0.1430  0.3781     85     no             obs   no 
## 
## Test for Residual Heterogeneity:
## QE(df = 83) = 312.3189, p-val < .0001
## 
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 83) = 4.2481, p-val = 0.0424
## 
## Model Results:
## 
##            estimate      se     tval  df    pval    ci.lb   ci.ub 
## intrcpt     -0.1036  0.6014  -0.1722  83  0.8637  -1.2997  1.0926    
## SSD_index    1.4134  0.6858   2.0611  83  0.0424   0.0495  2.7774  * 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    #lnCVR
    MLMR_mods_pers_SSD_mammal_aggression_lncvr <- rma.mv(CVR_yi ~ SSD_index, V = CVR_vi, 
                                            random = list(~1|study_ID, ~1|spp_name_phylo, ~1|obs), 
                                            R = list(spp_name_phylo=phylo_vcv_mammal), control=list(optimizer="optim"), 
                                            test = "t", data = pers_new_mammal_aggression)
## Warning in rma.mv(CVR_yi ~ SSD_index, V = CVR_vi, random = list(~1 | study_ID, :
## There are rows/columns in the 'R' matrix for 'spp_name_phylo' for which there
## are no data.

## Warning in rma.mv(CVR_yi ~ SSD_index, V = CVR_vi, random = list(~1 | study_ID, :
## Rows with NAs omitted from model fitting.
    MLMR_mods_pers_SSD_mammal_aggression_lncvr
## 
## Multivariate Meta-Analysis Model (k = 85; method: REML)
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed          factor    R 
## sigma^2.1  0.1790  0.4231     15     no        study_ID   no 
## sigma^2.2  0.0000  0.0052     13     no  spp_name_phylo  yes 
## sigma^2.3  0.1539  0.3922     85     no             obs   no 
## 
## Test for Residual Heterogeneity:
## QE(df = 83) = 202.3514, p-val < .0001
## 
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 83) = 0.0118, p-val = 0.9138
## 
## Model Results:
## 
##            estimate      se     tval  df    pval    ci.lb   ci.ub 
## intrcpt      0.0992  0.1534   0.6467  83  0.5196  -0.2059  0.4042    
## SSD_index   -0.0756  0.6959  -0.1086  83  0.9138  -1.4596  1.3084    
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Exploration:

  # d. exploration
     #SMD
    MLMR_mods_pers_SSD_mammal_explore_SMD <- rma.mv(SMD_yi_flip ~ SSD_index, V = SMD_vi, 
                                          random = list(~1|study_ID, ~1|spp_name_phylo, ~1|obs), 
                                          R = list(spp_name_phylo=phylo_vcv_mammal), control=list(optimizer="optim"), 
                                          test = "t", data = pers_new_mammal_exploration)
## Warning in rma.mv(SMD_yi_flip ~ SSD_index, V = SMD_vi, random = list(~1 | :
## There are rows/columns in the 'R' matrix for 'spp_name_phylo' for which there
## are no data.
    MLMR_mods_pers_SSD_mammal_explore_SMD
## 
## Multivariate Meta-Analysis Model (k = 213; method: REML)
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed          factor    R 
## sigma^2.1  0.0504  0.2244     19     no        study_ID   no 
## sigma^2.2  0.0000  0.0048     16     no  spp_name_phylo  yes 
## sigma^2.3  0.1331  0.3649    213     no             obs   no 
## 
## Test for Residual Heterogeneity:
## QE(df = 211) = 658.4587, p-val < .0001
## 
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 211) = 0.0351, p-val = 0.8516
## 
## Model Results:
## 
##            estimate      se     tval   df    pval    ci.lb   ci.ub 
## intrcpt     -0.0016  0.0914  -0.0173  211  0.9862  -0.1817  0.1786    
## SSD_index   -0.0522  0.2786  -0.1873  211  0.8516  -0.6015  0.4971    
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    #lnCVR
    MLMR_mods_pers_SSD_mammal_explore_lncvr <- rma.mv(CVR_yi ~ SSD_index, V = CVR_vi, 
                                            random = list(~1|study_ID, ~1|spp_name_phylo, ~1|obs), 
                                            R = list(spp_name_phylo=phylo_vcv_mammal), control=list(optimizer="optim"), 
                                            test = "t", data = pers_new_mammal_exploration)
## Warning in rma.mv(CVR_yi ~ SSD_index, V = CVR_vi, random = list(~1 | study_ID, :
## There are rows/columns in the 'R' matrix for 'spp_name_phylo' for which there
## are no data.
    MLMR_mods_pers_SSD_mammal_explore_lncvr
## 
## Multivariate Meta-Analysis Model (k = 213; method: REML)
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed          factor    R 
## sigma^2.1  0.0198  0.1406     19     no        study_ID   no 
## sigma^2.2  0.0265  0.1628     16     no  spp_name_phylo  yes 
## sigma^2.3  0.0323  0.1799    213     no             obs   no 
## 
## Test for Residual Heterogeneity:
## QE(df = 211) = 361.1620, p-val < .0001
## 
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 211) = 0.2658, p-val = 0.6067
## 
## Model Results:
## 
##            estimate      se     tval   df    pval    ci.lb   ci.ub 
## intrcpt     -0.0595  0.1507  -0.3951  211  0.6932  -0.3566  0.2375    
## SSD_index    0.1324  0.2567   0.5156  211  0.6067  -0.3737  0.6384    
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

8.5.3 Birds

SSD for boldness only.

# 2. BIRDS

  # subset dataset
  pers_new_bird <- as.data.frame(pers_new %>%
    filter(personality_trait == "boldness" & taxo_group == "bird"))

  # phylo_vcv birds only
  phylo_vcv_bird <- phylo_vcv[[1]]

Boldness:

    # SMD
    MLMR_mods_pers_SSD_bird_SMD <- rma.mv(SMD_yi_flip ~ SSD_index, V = SMD_vi, 
                                          random = list(~1|study_ID, ~1|spp_name_phylo, ~1|obs), 
                                          R = list(spp_name_phylo=phylo_vcv_bird), control=list(optimizer="optim"), 
                                          test = "t", data = pers_new_bird)
## Warning in rma.mv(SMD_yi_flip ~ SSD_index, V = SMD_vi, random = list(~1 | :
## There are rows/columns in the 'R' matrix for 'spp_name_phylo' for which there
## are no data.
## Warning: Rows with NAs omitted from model fitting.
    MLMR_mods_pers_SSD_bird_SMD
## 
## Multivariate Meta-Analysis Model (k = 234; method: REML)
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed          factor    R 
## sigma^2.1  1.9496  1.3963     21     no        study_ID   no 
## sigma^2.2  0.0001  0.0074     78     no  spp_name_phylo  yes 
## sigma^2.3  0.0925  0.3042    234     no             obs   no 
## 
## Test for Residual Heterogeneity:
## QE(df = 232) = 1579.6588, p-val < .0001
## 
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 232) = 0.1117, p-val = 0.7385
## 
## Model Results:
## 
##            estimate      se     tval   df    pval    ci.lb   ci.ub 
## intrcpt     -0.2211  0.3139  -0.7043  232  0.4819  -0.8397  0.3974    
## SSD_index   -0.2015  0.6028  -0.3342  232  0.7385  -1.3890  0.9861    
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    # lnCVR
    MLMR_mods_pers_SSD_bird_lncvr <- rma.mv(CVR_yi ~ SSD_index, V = CVR_vi, 
                                            random = list(~1|study_ID, ~1|spp_name_phylo, ~1|obs), 
                                            R = list(spp_name_phylo=phylo_vcv_bird), control=list(optimizer="optim"), 
                                            test = "t", data = pers_new_bird)
## Warning in rma.mv(CVR_yi ~ SSD_index, V = CVR_vi, random = list(~1 | study_ID, :
## There are rows/columns in the 'R' matrix for 'spp_name_phylo' for which there
## are no data.

## Warning in rma.mv(CVR_yi ~ SSD_index, V = CVR_vi, random = list(~1 | study_ID, :
## Rows with NAs omitted from model fitting.
    MLMR_mods_pers_SSD_bird_lncvr
## 
## Multivariate Meta-Analysis Model (k = 234; method: REML)
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed          factor    R 
## sigma^2.1  0.0000  0.0003     21     no        study_ID   no 
## sigma^2.2  0.0029  0.0537     78     no  spp_name_phylo  yes 
## sigma^2.3  0.0000  0.0003    234     no             obs   no 
## 
## Test for Residual Heterogeneity:
## QE(df = 232) = 244.9667, p-val = 0.2670
## 
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 232) = 0.6126, p-val = 0.4346
## 
## Model Results:
## 
##            estimate      se    tval   df    pval    ci.lb   ci.ub 
## intrcpt      0.0428  0.0361  1.1851  232  0.2372  -0.0283  0.1139    
## SSD_index    0.1023  0.1307  0.7827  232  0.4346  -0.1553  0.3599    
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

8.5.4 Fish

SSD for aggression and boldness.

# 3. FISH
  
  # subset by trait type
  # a. aggression
    pers_new_fish_aggression <- as.data.frame(pers_new %>%
    filter(personality_trait == "aggression") %>%
    filter(taxo_group == "fish")) 
  # b. boldness
    pers_new_fish_bold <- as.data.frame(pers_new %>%
    filter(personality_trait == "boldness") %>%
    filter(taxo_group == "fish")) 
    
  # phylo
    phylo_vcv_fish <- phylo_vcv[[2]]

Aggression:

  # a. aggression
    # SMD
      MLMR_mods_pers_SSD_fish_aggression_SMD <- rma.mv(SMD_yi_flip ~ SSD_index, V = SMD_vi, 
                                          random = list(~1|study_ID, ~1|spp_name_phylo, ~1|obs), 
                                          R = list(spp_name_phylo=phylo_vcv_fish), control=list(optimizer="optim"), 
                                          test = "t", data = pers_new_fish_aggression)
## Warning in rma.mv(SMD_yi_flip ~ SSD_index, V = SMD_vi, random = list(~1 | :
## There are rows/columns in the 'R' matrix for 'spp_name_phylo' for which there
## are no data.
## Warning: Rows with NAs omitted from model fitting.
    MLMR_mods_pers_SSD_fish_aggression_SMD
## 
## Multivariate Meta-Analysis Model (k = 93; method: REML)
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed          factor    R 
## sigma^2.1  0.0194  0.1395     16     no        study_ID   no 
## sigma^2.2  0.3329  0.5770     13     no  spp_name_phylo  yes 
## sigma^2.3  0.1704  0.4128     93     no             obs   no 
## 
## Test for Residual Heterogeneity:
## QE(df = 91) = 334.1728, p-val < .0001
## 
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 91) = 0.2301, p-val = 0.6326
## 
## Model Results:
## 
##            estimate      se     tval  df    pval    ci.lb   ci.ub 
## intrcpt     -0.1643  0.3987  -0.4120  91  0.6813  -0.9562  0.6277    
## SSD_index    0.2659  0.5544   0.4797  91  0.6326  -0.8352  1.3671    
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    # lnCVR
      MLMR_mods_pers_SSD_fish_aggression_lncvr <- rma.mv(CVR_yi ~ SSD_index, V = CVR_vi, 
                                            random = list(~1|study_ID, ~1|spp_name_phylo, ~1|obs), 
                                            R = list(spp_name_phylo=phylo_vcv_fish), control=list(optimizer="optim"), 
                                            test = "t", data = pers_new_fish_aggression)
## Warning in rma.mv(CVR_yi ~ SSD_index, V = CVR_vi, random = list(~1 | study_ID, :
## There are rows/columns in the 'R' matrix for 'spp_name_phylo' for which there
## are no data.

## Warning in rma.mv(CVR_yi ~ SSD_index, V = CVR_vi, random = list(~1 | study_ID, :
## Rows with NAs omitted from model fitting.
    MLMR_mods_pers_SSD_fish_aggression_lncvr
## 
## Multivariate Meta-Analysis Model (k = 93; method: REML)
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed          factor    R 
## sigma^2.1  0.0210  0.1450     16     no        study_ID   no 
## sigma^2.2  0.0000  0.0021     13     no  spp_name_phylo  yes 
## sigma^2.3  0.0000  0.0012     93     no             obs   no 
## 
## Test for Residual Heterogeneity:
## QE(df = 91) = 68.2701, p-val = 0.9640
## 
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 91) = 0.1495, p-val = 0.6999
## 
## Model Results:
## 
##            estimate      se     tval  df    pval    ci.lb   ci.ub 
## intrcpt     -0.1163  0.0597  -1.9483  91  0.0545  -0.2348  0.0023  . 
## SSD_index   -0.1323  0.3423  -0.3866  91  0.6999  -0.8122  0.5476    
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Boldness:

    # SMD
      MLMR_mods_pers_SSD_fish_bold_SMD <- rma.mv(SMD_yi_flip ~ SSD_index, V = SMD_vi, 
                                          random = list(~1|study_ID, ~1|spp_name_phylo, ~1|obs), 
                                          R = list(spp_name_phylo=phylo_vcv_fish), control=list(optimizer="optim"), 
                                          test = "t", data = pers_new_fish_bold)
## Warning in rma.mv(SMD_yi_flip ~ SSD_index, V = SMD_vi, random = list(~1 | :
## There are rows/columns in the 'R' matrix for 'spp_name_phylo' for which there
## are no data.
## Warning: Rows with NAs omitted from model fitting.
    MLMR_mods_pers_SSD_fish_bold_SMD
## 
## Multivariate Meta-Analysis Model (k = 171; method: REML)
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed          factor    R 
## sigma^2.1  0.1717  0.4143     23     no        study_ID   no 
## sigma^2.2  0.0279  0.1671     12     no  spp_name_phylo  yes 
## sigma^2.3  0.1634  0.4042    171     no             obs   no 
## 
## Test for Residual Heterogeneity:
## QE(df = 169) = 614.1157, p-val < .0001
## 
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 169) = 0.3196, p-val = 0.5726
## 
## Model Results:
## 
##            estimate      se     tval   df    pval    ci.lb   ci.ub 
## intrcpt      0.1883  0.1764   1.0678  169  0.2871  -0.1599  0.5366    
## SSD_index   -0.2571  0.4548  -0.5653  169  0.5726  -1.1550  0.6408    
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    # lnCVR
      MLMR_mods_pers_SSD_fish_bold_lncvr <- rma.mv(CVR_yi ~ SSD_index, V = CVR_vi, 
                                            random = list(~1|study_ID, ~1|spp_name_phylo, ~1|obs), 
                                            R = list(spp_name_phylo=phylo_vcv_fish), control=list(optimizer="optim"), 
                                            test = "t", data = pers_new_fish_bold)  
## Warning in rma.mv(CVR_yi ~ SSD_index, V = CVR_vi, random = list(~1 | study_ID, :
## There are rows/columns in the 'R' matrix for 'spp_name_phylo' for which there
## are no data.

## Warning in rma.mv(CVR_yi ~ SSD_index, V = CVR_vi, random = list(~1 | study_ID, :
## Rows with NAs omitted from model fitting.
    MLMR_mods_pers_SSD_fish_bold_lncvr
## 
## Multivariate Meta-Analysis Model (k = 171; method: REML)
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed          factor    R 
## sigma^2.1  0.0445  0.2109     23     no        study_ID   no 
## sigma^2.2  0.0184  0.1356     12     no  spp_name_phylo  yes 
## sigma^2.3  0.0881  0.2968    171     no             obs   no 
## 
## Test for Residual Heterogeneity:
## QE(df = 169) = 395.3375, p-val < .0001
## 
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 169) = 0.3068, p-val = 0.5804
## 
## Model Results:
## 
##            estimate      se     tval   df    pval    ci.lb   ci.ub 
## intrcpt      0.0063  0.1310   0.0482  169  0.9616  -0.2523  0.2649    
## SSD_index   -0.1507  0.2721  -0.5539  169  0.5804  -0.6878  0.3864    
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

8.5.5 Inverts

SSD for activity and boldness.

# 4. INVERTS
  # subset dataset
  # a. activity
    invert_activity <- as.data.frame(pers_new %>%
    filter(personality_trait == "activity") %>%
    filter(taxo_group == "invertebrate")) 
  # b. boldness
    invert_bold <- as.data.frame(pers_new %>%
    filter(personality_trait == "boldness") %>%
    filter(taxo_group == "invertebrate"))
  
  # phylo
    phylo_vcv_invert <- phylo_vcv[[3]]

Activity:

  # rerun models
    # a. activity
    # SMD
      MLMR_mods_pers_SSD_invert_activity_SMD <- rma.mv(SMD_yi_flip ~ SSD_index, V = SMD_vi, 
                                          random = list(~1|study_ID, ~1|spp_name_phylo, ~1|obs), 
                                          R = list(spp_name_phylo=phylo_vcv_invert), control=list(optimizer="optim"), 
                                          test = "t", data = invert_activity)
## Warning in rma.mv(SMD_yi_flip ~ SSD_index, V = SMD_vi, random = list(~1 | :
## There are rows/columns in the 'R' matrix for 'spp_name_phylo' for which there
## are no data.
## Warning: Rows with NAs omitted from model fitting.
    MLMR_mods_pers_SSD_invert_activity_SMD
## 
## Multivariate Meta-Analysis Model (k = 165; method: REML)
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed          factor    R 
## sigma^2.1  2.1874  1.4790     18     no        study_ID   no 
## sigma^2.2  0.0001  0.0104     16     no  spp_name_phylo  yes 
## sigma^2.3  0.1562  0.3952    165     no             obs   no 
## 
## Test for Residual Heterogeneity:
## QE(df = 163) = 1081.7241, p-val < .0001
## 
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 163) = 0.7140, p-val = 0.3993
## 
## Model Results:
## 
##            estimate      se     tval   df    pval    ci.lb   ci.ub 
## intrcpt      0.3479  0.3670   0.9480  163  0.3445  -0.3767  1.0725    
## SSD_index   -0.6862  0.8120  -0.8450  163  0.3993  -2.2896  0.9173    
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    # lnCVR
      MLMR_mods_pers_SSD_invert_activity_lncvr <- rma.mv(CVR_yi ~ SSD_index, V = CVR_vi, 
                                            random = list(~1|study_ID, ~1|spp_name_phylo, ~1|obs), 
                                            R = list(spp_name_phylo=phylo_vcv_invert), control=list(optimizer="optim"), 
                                            test = "t", data = invert_activity) 
## Warning in rma.mv(CVR_yi ~ SSD_index, V = CVR_vi, random = list(~1 | study_ID, :
## There are rows/columns in the 'R' matrix for 'spp_name_phylo' for which there
## are no data.

## Warning in rma.mv(CVR_yi ~ SSD_index, V = CVR_vi, random = list(~1 | study_ID, :
## Rows with NAs omitted from model fitting.
      MLMR_mods_pers_SSD_invert_activity_lncvr
## 
## Multivariate Meta-Analysis Model (k = 165; method: REML)
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed          factor    R 
## sigma^2.1  0.1193  0.3454     18     no        study_ID   no 
## sigma^2.2  0.0000  0.0035     16     no  spp_name_phylo  yes 
## sigma^2.3  0.0591  0.2431    165     no             obs   no 
## 
## Test for Residual Heterogeneity:
## QE(df = 163) = 486.7410, p-val < .0001
## 
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 163) = 0.4392, p-val = 0.5085
## 
## Model Results:
## 
##            estimate      se     tval   df    pval    ci.lb   ci.ub 
## intrcpt     -0.0278  0.1074  -0.2589  163  0.7960  -0.2400  0.1843    
## SSD_index    0.2685  0.4052   0.6627  163  0.5085  -0.5315  1.0685    
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Boldness:

   # SMD
      MLMR_mods_pers_SSD_invert_bold_SMD <- rma.mv(SMD_yi_flip ~ SSD_index, V = SMD_vi, 
                                          random = list(~1|study_ID, ~1|spp_name_phylo, ~1|obs), 
                                          R = list(spp_name_phylo=phylo_vcv_invert), control=list(optimizer="optim"), 
                                          test = "t", data = invert_bold)
## Warning in rma.mv(SMD_yi_flip ~ SSD_index, V = SMD_vi, random = list(~1 | :
## There are rows/columns in the 'R' matrix for 'spp_name_phylo' for which there
## are no data.
    MLMR_mods_pers_SSD_invert_bold_SMD
## 
## Multivariate Meta-Analysis Model (k = 164; method: REML)
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed          factor    R 
## sigma^2.1  0.0822  0.2867     23     no        study_ID   no 
## sigma^2.2  0.0000  0.0020     23     no  spp_name_phylo  yes 
## sigma^2.3  0.0650  0.2550    164     no             obs   no 
## 
## Test for Residual Heterogeneity:
## QE(df = 162) = 513.4222, p-val < .0001
## 
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 162) = 0.1533, p-val = 0.6959
## 
## Model Results:
## 
##            estimate      se    tval   df    pval    ci.lb   ci.ub 
## intrcpt      0.0985  0.0823  1.1967  162  0.2332  -0.0640  0.2611    
## SSD_index    0.1313  0.3354  0.3915  162  0.6959  -0.5310  0.7936    
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    # lnCVR
      MLMR_mods_pers_SSD_invert_bold_lncvr <- rma.mv(CVR_yi ~ SSD_index, V = CVR_vi, 
                                            random = list(~1|study_ID, ~1|spp_name_phylo, ~1|obs), 
                                            R = list(spp_name_phylo=phylo_vcv_invert), control=list(optimizer="optim"), 
                                            test = "t", data = invert_bold)  
## Warning in rma.mv(CVR_yi ~ SSD_index, V = CVR_vi, random = list(~1 | study_ID, :
## There are rows/columns in the 'R' matrix for 'spp_name_phylo' for which there
## are no data.
    MLMR_mods_pers_SSD_invert_bold_lncvr
## 
## Multivariate Meta-Analysis Model (k = 164; method: REML)
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed          factor    R 
## sigma^2.1  0.0383  0.1957     23     no        study_ID   no 
## sigma^2.2  0.0000  0.0015     23     no  spp_name_phylo  yes 
## sigma^2.3  0.0378  0.1945    164     no             obs   no 
## 
## Test for Residual Heterogeneity:
## QE(df = 162) = 380.7049, p-val < .0001
## 
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 162) = 0.0013, p-val = 0.9716
## 
## Model Results:
## 
##            estimate      se     tval   df    pval    ci.lb   ci.ub 
## intrcpt     -0.0145  0.0607  -0.2384  162  0.8119  -0.1344  0.1055    
## SSD_index   -0.0089  0.2503  -0.0356  162  0.9716  -0.5032  0.4854    
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

8.6 Multiple testing

Now we can extract the p-values from our intercept models, personality trait models, and SSD subset models to adjust p-values using the false discovery rate method. This method uses the p.adjust function to adjust p-values to account for multiple testing.

## Extract p-values from SSD subset models (reported in main text)

# list
  p.SMD_SSD <- c(0.69, 0.03, 0.86, 0.04, 0.34, 0.25, 0.99, 0.85, 0.48, 0.74, 0.68, 0.63, 0.29, 0.57, 0.34, 0.40, 0.23, 0.70)
  p.lnCVR_SSD <- c(0.60, 0.72, 0.52, 0.91, 0.82, 0.15, 0.69, 0.61, 0.24, 0.43, 0.05, 0.70, 0.96, 0.58, 0.80, 0.51, 0.81, 0.97)
  
  # p adjustment on our 3 hypothesis-testing models
  
    #SMD
  p.adjust(p = c(p.SMD_intercept, p.SMD_pers, p.SMD_SSD), method = "fdr") 
##          bird          fish  invertebrate        mammal      reptilia 
##     0.8533333     0.8533333     0.8533333     0.8708500     0.8533333 
##         bird1         bird2         bird3         bird4         bird5 
##     0.8708500     0.8533333     0.8533333     0.8708500     0.4800000 
##         fish1         fish2         fish3         fish4         fish5 
##     0.8708500     0.8533333     0.8708500     0.8533333     0.8708500 
## invertebrate1 invertebrate2 invertebrate3 invertebrate4 invertebrate5 
##     0.8533333     0.8533333     0.8533333     0.9195744     0.8533333 
##       mammal1       mammal2       mammal3       mammal4       mammal5 
##     0.8533333     0.8708500     0.8533333     0.9292872     0.8708500 
##     reptilia1     reptilia2     reptilia3     reptilia4     reptilia5 
##     0.9195744     0.8533333     0.8533333     0.4800000     0.9900000 
##                                                                       
##     0.8708500     0.4800000     0.9195744     0.4800000     0.8533333 
##                                                                       
##     0.8533333     0.9900000     0.9195744     0.8533333     0.8708500 
##                                                                       
##     0.8708500     0.8708500     0.8533333     0.8708500     0.8533333 
##                                           
##     0.8533333     0.8533333     0.8708500
    #lncvr
  p.adjust(p = c(p.lnCVR_intercept, p.lnCVR_pers, p.lnCVR_SSD), method = "fdr") 
##          bird          fish  invertebrate        mammal      reptilia 
##     0.9513857     0.9513857     0.9513857     0.9513857     0.9513857 
##         bird1         bird2         bird3         bird4         bird5 
##     0.9513857     0.9513857     0.9513857     0.4948879     0.9513857 
##         fish1         fish2         fish3         fish4         fish5 
##     0.9513857     0.8505594     0.9513857     0.9513857     0.4948879 
## invertebrate1 invertebrate2 invertebrate3 invertebrate4 invertebrate5 
##     0.9513857     0.9513857     0.9513857     0.9513857     0.9513857 
##       mammal1       mammal2       mammal3       mammal4       mammal5 
##     0.9513857     0.9513857     0.9513857     0.9513857     0.9513857 
##     reptilia1     reptilia2     reptilia3     reptilia4     reptilia5 
##     0.9513857     0.7414563     0.9513857     0.9513857     0.9513857 
##                                                                       
##     0.9513857     0.9513857     0.9513857     0.9513857     0.9513857 
##                                                                       
##     0.9513857     0.9513857     0.9513857     0.9513857     0.9513857 
##                                                                       
##     0.7414563     0.9513857     0.9700000     0.9513857     0.9513857 
##                                           
##     0.9513857     0.9513857     0.9700000
  # these p-values are in the order presented in tables, so easy to replace old p-values with new ones

9 Exploratory analyses

We collected some additional information from the literature (mating system) and from studies that we expected would influence sex differences. These analyses are strictly exploratory and just compare categorical moderator terms.

9.1 mating system

Do effect sizes from monogamous or multiply-mating species differ? Model summaries presented in Supplementary Table S3.

# what have we got to work with?
    pers_new %>%
    group_by(taxo_group, mating_system) %>%
    filter(!is.na(mating_system))%>%
    summarise(n = n(), studies = length(unique(study_ID)), species = length(unique(spp_name_phylo))) # make a table of numbers
## `summarise()` has grouped output by 'taxo_group'. You can override using the `.groups` argument.
# reload model output
rerun_models == FALSE
## [1] TRUE
    if(rerun_models == TRUE){
      MLMR_models_pers_mating_system <- meta_model_fits(pers_new, phylo_vcv, type = "pers_mate")
      saveRDS(MLMR_models_pers_mating_system, "./output/MLMR_models_pers_mating_system")
    } else{
     MLMR_models_pers_mating_system <- readRDS("./output/MLMR_models_pers_mating_system")
    }

# Extract the SMD and lnCVR results
    smd_mods_mating_system <- MLMR_models_pers_mating_system["SMD",]
      
    lnCVR_mods_mating_system <- MLMR_models_pers_mating_system["lnCVR",]

9.2 age

Do effect sizes from adults (sexually mature) or juveniles differ? Model summaries presented in Supplementary Table S4

# make a table
  data.frame(pers_new %>%
  group_by(taxo_group, age) %>%
  summarise(n= n(), N_spp = length(unique(spp_name_phylo)), N_studies = length(unique(study_ID))))
## `summarise()` has grouped output by 'taxo_group'. You can override using the `.groups` argument.
# reload model output:
rerun_models == FALSE
## [1] TRUE
  if(rerun_models == TRUE){
      MLMR_models_pers_age <- meta_model_fits(pers_new, phylo_vcv, type = "age")
      saveRDS(MLMR_models_pers_age, "./output/MLMR_models_pers_age")
    } else{
     MLMR_models_pers_age <- readRDS("./output/MLMR_models_pers_age")
    }


# Extract the SMD and lnCVR results
  smd_mods_pers_age <- MLMR_models_pers_age["SMD",]
    
  lnCVR_mods_pers_age <- MLMR_models_pers_age["lnCVR",]

9.3 population

Do effect sizes from wild animals or lab animals differ? Model summaries presented in Supplementary Table S5.

# table
  data.frame(pers_new %>%
  group_by(taxo_group, population) %>%
  summarise(n = n(), N_spp = length(unique(spp_name_phylo)), N_studies = length(unique(study_ID))))
## `summarise()` has grouped output by 'taxo_group'. You can override using the `.groups` argument.
# reload model output:
rerun_models == FALSE
## [1] TRUE
  if(rerun_models == TRUE){
      MLMR_models_pers_pop <- meta_model_fits(pers_new, phylo_vcv, type = "pop")
      saveRDS(MLMR_models_pers_pop, "./output/MLMR_models_pers_pop")
    } else{
     MLMR_models_pers_pop <- readRDS("./output/MLMR_models_pers_pop")
    }

# Extract the SMD and lnCVR results
  smd_mods_pers_pop <- MLMR_models_pers_pop["SMD",]
    
  lnCVR_mods_pers_pop <- MLMR_models_pers_pop["lnCVR",]

9.4 study environment

Do effect sizes collected in the wild or the lab differ? Model summaries presented in Supplementary Table S6.

# table
  data.frame(pers_new %>%
  group_by(taxo_group, study_environment) %>%
  summarise(n = n(), N_spp = length(unique(spp_name_phylo)), N_studies = length(unique(study_ID))))
## `summarise()` has grouped output by 'taxo_group'. You can override using the `.groups` argument.
# reload model output:
rerun_models == FALSE
## [1] TRUE
  if(rerun_models == TRUE){
      MLMR_models_pers_environ <- meta_model_fits(pers_new, phylo_vcv, type = "environ")
      saveRDS(MLMR_models_pers_environ, "./output/MLMR_models_pers_environ")
    } else{
     MLMR_models_pers_environ <- readRDS("./output/MLMR_models_pers_environ")
    }

# Extract the SMD and lnCVR results
  smd_mods_pers_enviro <- MLMR_models_pers_environ["SMD",]
      
  lnCVR_mods_pers_enviro <- MLMR_models_pers_environ["lnCVR",]

9.5 study type

Do effect sizes from observational or experimental study design differ? Model summaries presented in Supplementary Table S7.

# Let's see what we have to work with
  data.frame(pers_new %>%
  group_by(taxo_group, study_type) %>%
  summarise(N_spp = length(unique(spp_name_phylo)), N_studies = length(unique(study_ID))))
## `summarise()` has grouped output by 'taxo_group'. You can override using the `.groups` argument.
    # inverts only have experimental observations, so need to exclude inverts from this analysis
    # because our phylo_vcv matrix is in a list that is hard to drop elements from, let's just run each model individually 

# 1. Mammals
  # Subset data
    pers_new_mammal <- as.data.frame(pers_new %>%
    filter(taxo_group == "mammal"))

  # Run models with just study type as moderator:
    #SMD
    MLMR_mods_pers_studytype_mammal_SMD <- rma.mv(SMD_yi_flip ~ study_type, V = SMD_vi, 
                                                   random = list(~1|study_ID, ~1|spp_name_phylo, ~1|obs), 
                                                   R = list(spp_name_phylo=phylo_vcv_mammal), control=list(optimizer="optim"), 
                                                   test = "t", data = pers_new_mammal)

    MLMR_mods_pers_studytype_mammal_SMD
## 
## Multivariate Meta-Analysis Model (k = 674; method: REML)
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed          factor    R 
## sigma^2.1  0.1051  0.3242     61     no        study_ID   no 
## sigma^2.2  0.0094  0.0971     45     no  spp_name_phylo  yes 
## sigma^2.3  0.1570  0.3963    674     no             obs   no 
## 
## Test for Residual Heterogeneity:
## QE(df = 672) = 2198.5222, p-val < .0001
## 
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 672) = 9.6851, p-val = 0.0019
## 
## Model Results:
## 
##                        estimate      se     tval   df    pval    ci.lb   ci.ub 
## intrcpt                 -0.0067  0.0913  -0.0734  672  0.9415  -0.1860  0.1726 
## study_typeobservation    0.4092  0.1315   3.1121  672  0.0019   0.1510  0.6674 
##  
## intrcpt 
## study_typeobservation  ** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
    #lnCVR
    MLMR_mods_pers_studytype_mammal_lncvr <- rma.mv(CVR_yi ~ study_type, V = CVR_vi, 
                                                   random = list(~1|study_ID, ~1|spp_name_phylo, ~1|obs), 
                                                   R = list(spp_name_phylo=phylo_vcv_mammal), control=list(optimizer="optim"), 
                                                   test = "t", data = pers_new_mammal)
    
    MLMR_mods_pers_studytype_mammal_lncvr
## 
## Multivariate Meta-Analysis Model (k = 674; method: REML)
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed          factor    R 
## sigma^2.1  0.0356  0.1888     61     no        study_ID   no 
## sigma^2.2  0.0436  0.2088     45     no  spp_name_phylo  yes 
## sigma^2.3  0.0338  0.1837    674     no             obs   no 
## 
## Test for Residual Heterogeneity:
## QE(df = 672) = 1061.6294, p-val < .0001
## 
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 672) = 0.5661, p-val = 0.4521
## 
## Model Results:
## 
##                        estimate      se    tval   df    pval    ci.lb   ci.ub 
## intrcpt                  0.0298  0.1361  0.2188  672  0.8269  -0.2375  0.2970 
## study_typeobservation    0.0771  0.1024  0.7524  672  0.4521  -0.1240  0.2781 
##  
## intrcpt 
## study_typeobservation 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# 2. BIRDS
  # subset dataset
      pers_new_bird <- as.data.frame(pers_new %>%
      filter(taxo_group == "bird"))

  # rerun models
      #SMD
       MLMR_mods_pers_studytype_bird_SMD <- rma.mv(SMD_yi_flip ~ study_type, V = SMD_vi, 
                                                   random = list(~1|study_ID, ~1|spp_name_phylo, ~1|obs), 
                                                   R = list(spp_name_phylo=phylo_vcv_bird), control=list(optimizer="optim"), 
                                                   test = "t", data = pers_new_bird)

       MLMR_mods_pers_studytype_bird_SMD
## 
## Multivariate Meta-Analysis Model (k = 480; method: REML)
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed          factor    R 
## sigma^2.1  0.6650  0.8155     50     no        study_ID   no 
## sigma^2.2  0.0000  0.0038    106     no  spp_name_phylo  yes 
## sigma^2.3  0.1203  0.3469    480     no             obs   no 
## 
## Test for Residual Heterogeneity:
## QE(df = 478) = 2378.8387, p-val < .0001
## 
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 478) = 1.0068, p-val = 0.3162
## 
## Model Results:
## 
##                        estimate      se     tval   df    pval    ci.lb   ci.ub 
## intrcpt                 -0.0210  0.1515  -0.1389  478  0.8896  -0.3187  0.2766 
## study_typeobservation   -0.2540  0.2532  -1.0034  478  0.3162  -0.7515  0.2434 
##  
## intrcpt 
## study_typeobservation 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
      #lnCVR
      MLMR_mods_pers_studytype_bird_lncvr <- rma.mv(CVR_yi ~ study_type, V = CVR_vi, 
                                                   random = list(~1|study_ID, ~1|spp_name_phylo, ~1|obs), 
                                                   R = list(spp_name_phylo=phylo_vcv_bird), control=list(optimizer="optim"), 
                                                   test = "t", data = pers_new_bird)

      MLMR_mods_pers_studytype_bird_lncvr
## 
## Multivariate Meta-Analysis Model (k = 480; method: REML)
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed          factor    R 
## sigma^2.1  0.2537  0.5036     50     no        study_ID   no 
## sigma^2.2  0.0001  0.0086    106     no  spp_name_phylo  yes 
## sigma^2.3  0.3766  0.6137    480     no             obs   no 
## 
## Test for Residual Heterogeneity:
## QE(df = 478) = 3186.9869, p-val < .0001
## 
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 478) = 1.7076, p-val = 0.1919
## 
## Model Results:
## 
##                        estimate      se     tval   df    pval    ci.lb   ci.ub 
## intrcpt                  0.0436  0.1109   0.3929  478  0.6946  -0.1744  0.2615 
## study_typeobservation   -0.2451  0.1876  -1.3067  478  0.1919  -0.6138  0.1235 
##  
## intrcpt 
## study_typeobservation 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# 3. FISH
  # subset dataset
      pers_new_fish <- as.data.frame(pers_new %>%
      filter(taxo_group == "fish")) 

  # rerun models
      #SMD
      MLMR_mods_pers_studytype_fish_SMD <- rma.mv(SMD_yi_flip ~ study_type, V = SMD_vi, 
                                          random = list(~1|study_ID, ~1|spp_name_phylo, ~1|obs), 
                                          R = list(spp_name_phylo=phylo_vcv_fish), control=list(optimizer="optim"), 
                                          test = "t", data = pers_new_fish)

      MLMR_mods_pers_studytype_fish_SMD
## 
## Multivariate Meta-Analysis Model (k = 490; method: REML)
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed          factor    R 
## sigma^2.1  0.5997  0.7744     44     no        study_ID   no 
## sigma^2.2  0.0440  0.2098     22     no  spp_name_phylo  yes 
## sigma^2.3  0.1091  0.3303    490     no             obs   no 
## 
## Test for Residual Heterogeneity:
## QE(df = 488) = 1523.5807, p-val < .0001
## 
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 488) = 0.0188, p-val = 0.8911
## 
## Model Results:
## 
##                        estimate      se     tval   df    pval    ci.lb   ci.ub 
## intrcpt                  0.1825  0.2106   0.8669  488  0.3864  -0.2312  0.5962 
## study_typeobservation   -0.0657  0.4798  -0.1369  488  0.8911  -1.0085  0.8771 
##  
## intrcpt 
## study_typeobservation 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
      #lnCVR
      MLMR_mods_pers_studytype_fish_lncvr <- rma.mv(CVR_yi ~ study_type, V = CVR_vi, 
                                            random = list(~1|study_ID, ~1|spp_name_phylo, ~1|obs), 
                                            R = list(spp_name_phylo=phylo_vcv_fish), control=list(optimizer="optim"), 
                                            test = "t", data = pers_new_fish)

      MLMR_mods_pers_studytype_fish_lncvr
## 
## Multivariate Meta-Analysis Model (k = 490; method: REML)
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed          factor    R 
## sigma^2.1  0.0352  0.1875     44     no        study_ID   no 
## sigma^2.2  0.0017  0.0408     22     no  spp_name_phylo  yes 
## sigma^2.3  0.1094  0.3307    490     no             obs   no 
## 
## Test for Residual Heterogeneity:
## QE(df = 488) = 1122.6615, p-val < .0001
## 
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 488) = 0.0722, p-val = 0.7882
## 
## Model Results:
## 
##                        estimate      se     tval   df    pval    ci.lb   ci.ub 
## intrcpt                 -0.0025  0.0536  -0.0460  488  0.9633  -0.1078  0.1029 
## study_typeobservation   -0.0399  0.1486  -0.2687  488  0.7882  -0.3319  0.2520 
##  
## intrcpt 
## study_typeobservation 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# 4. Reptiles
  # subset dataset
      pers_new_reptile <- as.data.frame(pers_new %>%
      filter(taxo_group == "reptilia")) 
 
  # phylo
      phylo_vcv_reptile <- phylo_vcv[[5]]
  
  # rerun models
      #SMD
      MLMR_mods_pers_studytype_rep_SMD <- rma.mv(SMD_yi_flip ~ study_type, V = SMD_vi, 
                                                random = list(~1|study_ID, ~1|spp_name_phylo, ~1|obs), 
                                          R = list(spp_name_phylo=phylo_vcv_reptile), control=list(optimizer="optim"), 
                                          test = "t", data = pers_new_reptile)

      MLMR_mods_pers_studytype_rep_SMD      
## 
## Multivariate Meta-Analysis Model (k = 95; method: REML)
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed          factor    R 
## sigma^2.1  0.0000  0.0008     11     no        study_ID   no 
## sigma^2.2  0.0730  0.2702     10     no  spp_name_phylo  yes 
## sigma^2.3  0.0426  0.2063     95     no             obs   no 
## 
## Test for Residual Heterogeneity:
## QE(df = 93) = 159.2776, p-val < .0001
## 
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 93) = 3.4462, p-val = 0.0666
## 
## Model Results:
## 
##                        estimate      se     tval  df    pval    ci.lb   ci.ub 
## intrcpt                  0.1334  0.1531   0.8715  93  0.3857  -0.1706  0.4374 
## study_typeobservation   -0.5085  0.2739  -1.8564  93  0.0666  -1.0524  0.0354 
##  
## intrcpt 
## study_typeobservation  . 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
      #lnCVR
      MLMR_mods_pers_studytype_rep_lncvr <- rma.mv(CVR_yi ~ study_type, V = CVR_vi, 
                                            random = list(~1|study_ID, ~1|spp_name_phylo, ~1|obs), 
                                            R = list(spp_name_phylo=phylo_vcv_reptile), control=list(optimizer="optim"), 
                                            test = "t", data = pers_new_reptile)
      
      MLMR_mods_pers_studytype_rep_lncvr
## 
## Multivariate Meta-Analysis Model (k = 95; method: REML)
## 
## Variance Components:
## 
##             estim    sqrt  nlvls  fixed          factor    R 
## sigma^2.1  0.0000  0.0003     11     no        study_ID   no 
## sigma^2.2  0.0000  0.0012     10     no  spp_name_phylo  yes 
## sigma^2.3  0.0000  0.0003     95     no             obs   no 
## 
## Test for Residual Heterogeneity:
## QE(df = 93) = 58.0505, p-val = 0.9983
## 
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 93) = 0.1931, p-val = 0.6613
## 
## Model Results:
## 
##                        estimate      se     tval  df    pval    ci.lb   ci.ub 
## intrcpt                  0.0427  0.0420   1.0173  93  0.3116  -0.0407  0.1261 
## study_typeobservation   -0.0621  0.1413  -0.4395  93  0.6613  -0.3426  0.2184 
##  
## intrcpt 
## study_typeobservation 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

10 Sensitivity analyses - Dependency matrix models

We need to refit our 3 main models accounting for any dependency resulting from the same traits measured on the same animals (likely a big source of non-independence) and any other shared covariance. We added the D matrices to the residual variance matrix as opposed to the sampling covariance. We chose to set 3 different levels of dependency (rho): 0.3, 0.5 an 0.8.

Model summaries are also presented in Supplementary Tables S8-S13.

 # Create the dependency matrices; try 3 levels of rho = 0.3, 0.5, 0.8
      
      pers_new <- data.frame(pers_new %>%
      group_by(taxo_group) %>%
      mutate(depend_n = paste0(study_ID, "_", depend)))
    
      split_taxa <- split(pers_new, pers_new$taxo_group)
    
  # 0.3 rho:
    D_matrices_0.3 <- lapply(split_taxa, function(x) make_VCV_matrix(x, V = x$SMD_vi, cluster = "depend_n", 
                                                                     obs = "obs", type = "cor", rho = 0.3))
  # 0.5 rho:
    D_matrices_0.5 <- lapply(split_taxa, function(x) make_VCV_matrix(x, V = x$SMD_vi, cluster = "depend_n", 
                                                                     obs = "obs", type = "cor", rho = 0.5))
  # 0.8 rho:   
    D_matrices_0.8 <- lapply(split_taxa, function(x) make_VCV_matrix(x, V = x$SMD_vi, cluster = "depend_n", 
                                                                     obs = "obs", type = "cor", rho = 0.8))

10.1 Intercept-only models with D matrices

Model output is presented in Supplementary Tables S8-S10 in the Supporting Information.

# 1. Intercept only models
    # rho = 0.3
    int_0.3 <- fit_int_MLMAmodel_D(pers_new, phylo_vcv, D_matrices_0.3)
    
    smd_mods_D_0.3 <- int_0.3[["SMD"]] 
    lnCVR_mods_D_0.3 <- int_0.3[["lnCVR"]] 
  
   # prediction intervals
    MLMA_estimates_SMD_D_0.3 <- plyr::ldply(lapply(smd_mods_D_0.3, 
                           function(x) print(mod_results(x, mod = "Int")))) 
    MLMA_estimates_lnCVR_D_0.3 <- plyr::ldply(lapply(lnCVR_mods_D_0.3, 
                              function(x) print(mod_results(x, mod = "Int"))))   

    # rho = 0.5
    int_0.5 <- fit_int_MLMAmodel_D(pers_new, phylo_vcv, D_matrices_0.5)
    
    smd_mods_D_0.5 <- int_0.5[["SMD"]] 
    lnCVR_mods_D_0.5 <- int_0.5[["lnCVR"]] 
  
    # prediction intervals
    MLMA_estimates_SMD_D_0.5 <- plyr::ldply(lapply(smd_mods_D_0.5, 
                             function(x) print(mod_results(x, mod = "Int")))) 
    MLMA_estimates_lnCVR_D_0.5 <- plyr::ldply(lapply(lnCVR_mods_D_0.5, 
                              function(x) print(mod_results(x, mod = "Int"))))
    
    # rho = 0.8
    int_0.8 <- fit_int_MLMAmodel_D(pers_new, phylo_vcv, D_matrices_0.8) 
    
    smd_mods_D_0.8 <- int_0.8[["SMD"]] 
    lnCVR_mods_D_0.8 <- int_0.8[["lnCVR"]] 
    
    # prediction intervals
    MLMA_estimates_SMD_D_0.8 <- plyr::ldply(lapply(smd_mods_D_0.8, 
                           function(x) print(mod_results(x, mod = "Int")))) 
    MLMA_estimates_lnCVR_D_0.8 <- plyr::ldply(lapply(lnCVR_mods_D_0.8, 
                              function(x) print(mod_results(x, mod = "Int"))))

10.2 Personality trait models with D matrices

Model output is presented in Supplementary Tables S11-S13 in the Supporting Information.

# 2. Personality Trait models

  # rho = 0.3
    pers_0.3 <- fit_int_MLMAmodel_D_pers(pers_new, phylo_vcv, D_matrices_0.3)
    
    smd_mods_D_pers_0.3 <- pers_0.3[["SMD"]] 
    lnCVR_mods_D_pers_0.3 <- pers_0.3[["lnCVR"]] 
  
    # prediction intervals
    MLMA_estimates_SMD_pers_D_0.3 <- plyr::ldply(lapply(smd_mods_D_pers_0.3, 
                           function(x) print(mod_results(x, mod = "personality_trait")))) 
    MLMA_estimates_lnCVR_pers_D_0.3 <- plyr::ldply(lapply(lnCVR_mods_D_pers_0.3, 
                              function(x) print(mod_results(x, mod = "personality_trait")))) 

  # rho = 0.5
    pers_0.5 <- fit_int_MLMAmodel_D_pers(pers_new, phylo_vcv, D_matrices_0.5)
    
    smd_mods_D_pers_0.5 <- pers_0.5[["SMD"]] 
    lnCVR_mods_D_pers_0.5 <- pers_0.5[["lnCVR"]] 
  
    # prediction intervals
    MLMA_estimates_SMD_pers_D_0.5 <- plyr::ldply(lapply(smd_mods_D_pers_0.5, 
                           function(x) print(mod_results(x, mod = "personality_trait")))) 
    MLMA_estimates_lnCVR_pers_D_0.5 <- plyr::ldply(lapply(lnCVR_mods_D_pers_0.5, 
                              function(x) print(mod_results(x, mod = "personality_trait")))) 

  # rho = 0.8
    pers_0.8 <- fit_int_MLMAmodel_D_pers(pers_new, phylo_vcv, D_matrices_0.8)
    
    smd_mods_D_pers_0.8 <- pers_0.8[["SMD"]] 
    lnCVR_mods_D_pers_0.8 <- pers_0.8[["lnCVR"]] 
  
    # prediction intervals
    MLMA_estimates_SMD_pers_D_0.8 <- plyr::ldply(lapply(smd_mods_D_pers_0.8, 
                           function(x) print(mod_results(x, mod = "personality_trait")))) 
    MLMA_estimates_lnCVR_pers_D_0.8 <- plyr::ldply(lapply(lnCVR_mods_D_pers_0.8, 
                              function(x) print(mod_results(x, mod = "personality_trait")))) 

10.3 Personality * SSD models

These models were just to check since we don't really interpret the interaction models.

# 3. Pers Trait * SSD models 
# just use the full interaction models here since this is just a check
# won't bother with prediction intervals here since these models aren't really for interpretation

    # rho = 0.3
      ssd_0.3 <- fit_int_MLMAmodel_D_pers_ssd(pers_new, phylo_vcv, D_matrices_0.3)
      
      split_taxa <- split(pers_new, pers_new$taxo_group)
      smd_mods_D_pers_ssd_0.3 <- ssd_0.3[["SMD"]] 
      lnCVR_mods_D_pers_ssd_0.3 <- ssd_0.3[["lnCVR"]] 
  
    # rho = 0.5
      ssd_0.5 <- fit_int_MLMAmodel_D_pers_ssd(pers_new, phylo_vcv, D_matrices_0.5)
    
      split_taxa <- split(pers_new, pers_new$taxo_group)
      smd_mods_D_pers_ssd_0.5 <- ssd_0.5[["SMD"]] 
      lnCVR_mods_D_pers_ssd_0.5 <- ssd_0.5[["lnCVR"]]   
     
    # rho = 0.8
      ssd_0.8 <- fit_int_MLMAmodel_D_pers_ssd(pers_new, phylo_vcv, D_matrices_0.8)
    
      split_taxa <- split(pers_new, pers_new$taxo_group)
      smd_mods_D_pers_ssd_0.8 <- ssd_0.8[["SMD"]] 
      lnCVR_mods_D_pers_ssd_0.8 <- ssd_0.8[["lnCVR"]] 

11 Sensitivity analyses - Publication Bias

We can use: 1) funnel plots to look for asymmetry across all effect sizes for both SMD and lnCVR, and 2) Egger's test which performs a regression test on our funnel plots ... but is not useful when there is high heterogeneity NOT caused by publication bias (which is the case for our data).

Since our data has very high heterogeneity, we instead included the inverse of the 'effective sample size' as a moderator term in our full model (personality trait model) to see if study precision is driving effect size patterns. The logic here is that studies with low or high precision can have a significant influence and so including precision as a moderator will allow us to see if precision is significant (and which direction). See Nakagawa et al. 2021 for more info (reference in main text).

Model summaries are presented in Supplementary Table S18.

### NEW METHOD OF PUBLICATION BIAS FROM NAKAGAWA ET AL 2021 - PREPRINT

    # calculating the inverse of the 'effective sample size' to account for unbalanced sampling
    pers_new$inv_n_tilda <- with(pers_new, ((female_n + male_n)/(female_n*male_n)))
    pers_new$sqrt_inv_n_tilda <- with(pers_new, (sqrt(inv_n_tilda))) # use this in the model

    if(rerun_models == TRUE){
      MLMR_models_pers_pubbias <- meta_model_fits(pers_new, phylo_vcv, type = "pubbias")
      saveRDS(MLMR_models_pers_pubbias, "./output/MLMR_models_pers_pubbias")
    } else{
     MLMR_models_pers_pubbias <- readRDS("./output/MLMR_models_pers_pubbias")
    }

    # Extract the SMD and lnCVR results
    smd_mods_pubbias <- MLMR_models_pers_pubbias["SMD",] 
   
    lnCVR_mods_pubbias <- MLMR_models_pers_pubbias["lnCVR",] 

12 Exploratory analyses - Heterogamety and taxo group

There was a trend for male mammals to be more variable than females and female birds to be more variable than males, for some of the five personality traits. To better compare the direction of these effect sizes we decided post hoc to conduct an exploratory analysis with personality trait type and taxonomic group as moderator terms to compare birds and mammals (males homogametic or heterogametic, respectively). To do this, we first combined the bird and mammal phylo correlation matrices together (assuming no phylo heritability across the taxo groups - since phylo did not really explain heterogeneity it shouldn't matter). We then created an interaction MLMR model with personality trait * taxa (no intercept) to get slope estimates for each of the traits for mammals and birds seperately.

From this model, we then compared each of the five traits for birds and mammals using a post hoc Tukey pairwise comparison to test whether birds and mammals were significantly different from each other.

Model summaries are presented in Supplementary Table S19.

# install packages to make diagonal matrix and to make multiple comparisons 
  library(multcomp)
  library(Matrix)

# Create block diag phylogeny
  phylogeny <- Matrix::bdiag(phylo_vcv_bird, phylo_vcv_mammal) # use this as the phylo vcv in the model

    # needs to have colnames for use in random effects model
    dimnames(phylogeny) <- Map(c, dimnames(phylo_vcv_bird), dimnames(phylo_vcv_mammal))

    # only include bird and mammal data
    pers_new_contrast <- as.data.frame(pers_new %>%
                                        filter(taxo_group =="mammal" | taxo_group == "bird") %>% 
                                        mutate(sp_pers = interaction(personality_trait,taxo_group)))
    
    # 1. intercept only model
    contrast_birdmammal_lncvr_int <- rma.mv(CVR_yi ~ taxo_group, V = CVR_vi, 
                                                 random = list(~1|study_ID, ~1|spp_name_phylo, ~1|obs), 
                                                 R = list(spp_name_phylo=phylogeny), control=list(optimizer="optim"), 
                                                 test = "t", data = pers_new_contrast) 
    

    # 2. personality trait model  
    # creating the model - with pers trait and taxo group as mods 

    #lnCVR model only
    contrast_birdmammal_lncvr <- rma.mv(CVR_yi ~ personality_trait*taxo_group, V = CVR_vi, 
                                                 random = list(~1|study_ID, ~1|spp_name_phylo, ~1|obs), 
                                                 R = list(spp_name_phylo=phylogeny), control=list(optimizer="optim"), 
                                                 test = "t", data = pers_new_contrast)
    
    # model with interaction only to check output of model above
    contrast_birdmammal_lncvr_2 <- rma.mv(CVR_yi ~ sp_pers -1, V = CVR_vi, 
                                                 random = list(~1|study_ID, ~1|spp_name_phylo, ~1|obs), 
                                                 R = list(spp_name_phylo=phylogeny), control=list(optimizer="optim"), 
                                                 test = "t", data = pers_new_contrast)

# multiple comparison using Tukey test
  summary(glht(contrast_birdmammal_lncvr, linfct = cbind(contrMat(rep(1:10), type = "Tukey"))), test=adjusted("fdr"))

# here we are only interested in the comparisons between mammals and birds, so: 1-6 (activity), 2-7 (aggression), 3-8 (boldness), 4-9 (exploration), and 5-10 (sociality)

13 Plots

13.1 Orchard plots of effect sizes from personality trait models

These plots use the orchaRd package to generate pretty plots where each effect size (k) is a point on the plot.

13.1.1 lnCVR

# create objects of each of the models first

# lnCVR

  # Bird lnCVR
  bird_lncvr <- orchard_plot(lnCVR_mods_pers[[1]], mod = "personality_trait", xlab = "log Coefficient of Variance (lnCVR)", angle = 45, alpha = 0.5, transfm = "none")
  # Fish lnCVR
  fish_lncvr <- orchard_plot(lnCVR_mods_pers[[2]], mod = "personality_trait", xlab = "log Coefficient of Variance (lnCVR)", angle = 45, alpha = 0.5, transfm = "none")
  # Invert lnCVR
  invert_lncvr<- orchard_plot(lnCVR_mods_pers[[3]], mod = "personality_trait", xlab = "log Coefficient of Variance (lnCVR)", angle = 45, alpha = 0.5, transfm = "none")
  # Mammal lnCVR
  mammal_lncvr <- orchard_plot(lnCVR_mods_pers[[4]], mod = "personality_trait", xlab = "log Coefficient of Variance (lnCVR)", angle = 45, alpha = 0.5, transfm = "none")
  # Reptile lnCVR
  reptile_lncvr <- orchard_plot(lnCVR_mods_pers[[5]], mod = "personality_trait", xlab = "log Coefficient of Variance (lnCVR)", angle = 45, alpha = 0.5, transfm = "none")

13.1.2 SMD

  # Bird SMD
  bird_SMD <- orchard_plot(smd_mods_pers[[1]], mod = "personality_trait", xlab = "Standardised mean difference", angle = 45, alpha = 0.5, transfm = "none") 
  # Fish SMD
  fish_SMD <- orchard_plot(smd_mods_pers[[2]], mod = "personality_trait", xlab = "Standardised mean difference", angle = 45, alpha = 0.5, transfm = "none")
  # Invert SMD
  invert_SMD<- orchard_plot(smd_mods_pers[[3]], mod = "personality_trait", xlab = "Standardised mean difference", angle = 45, alpha = 0.5, transfm = "none")
  # Mammal SMD
  mammal_SMD <- orchard_plot(smd_mods_pers[[4]], mod = "personality_trait", xlab = "Standardised mean difference", angle = 45, alpha = 0.5, transfm = "none")
  # Reptile SMD
  reptile_SMD <- orchard_plot(smd_mods_pers[[5]], mod = "personality_trait", xlab = "Standardised mean difference", angle = 45, alpha = 0.5, transfm = "none")

Putting the SMD and lnCVR plots together

Endotherms:

  # window size for orchard plots
  # the precision guides on the plots are a bit ugly, collect them to the side and crop them out   
  
## Mammals
 
  # dev.new(width=8,height=7,noRStudioGD = TRUE)

# place plots together

mammal_SMD / mammal_lncvr 

# ggsave("./figs/finished figs/mammal_effects.tiff", width = 8, height = 7, units = "in") #save image

## Birds

bird_SMD / bird_lncvr

# ggsave("./figs/finished figs/bird_effects.tiff", width = 8, height = 7, units = "in") #save image

Ectotherms:

# window size a bit smaller for these guys
    # dev.new(width=8,height=5,noRStudioGD = TRUE)

## Reptiles and amphibians

reptile_SMD / reptile_lncvr / plot_layout(guides = 'collect')

# ggsave("~/Documents/GitHub/sex_meta/figs/finished figs/rep_effects.tiff", width = 8, height = 5, units = "in")

## Fish

fish_SMD / fish_lncvr / plot_layout(guides = 'collect')

# ggsave("~/Documents/GitHub/sex_meta/figs/finished figs/fish_effects.tiff", width = 8, height = 5, units = "in")

## Invertebrates

invert_SMD / invert_lncvr / plot_layout(guides = 'collect')

# ggsave("~/Documents/GitHub/sex_meta/figs/finished figs/invert_effects.tiff", width = 8, height = 5, units = "in")

The precision guides will get cropped out when joining the orchard plots and phylogenies together.


13.2 Phylogenetic trees with heatmaps

Using ggtree to plot lots of complex data onto phylogenetic trees see: https://guangchuangyu.github.io/ggtree-book/chapter-ggtree.html for more information about using ggtree

# install ggtree using this method:
  source("https://bioconductor.org/biocLite.R")
  BiocManager::install("ggtree")
## 
## The downloaded binary packages are in
##  /var/folders/0b/pxghylq157gfhs1vrzdpx2gc0000gq/T//RtmpXG0rME/downloaded_packages
  library(ggtree)

# load organised SSD data using figs_data.csv 
  figs_data <- read.csv("./data/figs_data.csv", stringsAsFactors = FALSE)

13.3 bird tree

# subset dataset to include only birds
 bird_data <- as.data.frame(figs_data %>%
    filter(taxo_group == "bird"))

# setting up the basic tree structure
  # load tree
    birdtree <- read.tree("./trees/bird_species.nwk")

  # prune tree to get rid of species we no longer have data for
    pruned.birdtree <- drop.tip(birdtree, setdiff(birdtree$tip.label, bird_data$spp_name_phylo)) 

  # remove underscores from tip labels
    pruned.birdtree$tip.label = gsub("_", " ", pruned.birdtree$tip.label)
   
  # remove underscores from species name in our species data list
    bird_data$spp_name_phylo = gsub("_", " ", bird_data$spp_name_phylo) 

  # set row names
    row.names(bird_data) <- bird_data$spp_name_phylo 

  # define objects for the plot
    species <- pruned.birdtree$tip.label

    rownames(bird_data) <- pruned.birdtree$tip.label 

# set window size
# dev.new(width=8, height=8,noRStudioGD = TRUE) #opens quartz window of set size

# now need to make a matrix of effect sizes (n) for each species for each personality trait to add to our plot!
    # subset dataset
     pers_bird <- as.data.frame(pers_new %>%
     filter(taxo_group == "bird")) 
  
    # make this a matrix-style dataframe
      pers_bird <- data.frame(pers_bird %>%
      group_by(spp_name_phylo, personality_trait) %>%
      summarise(n = n()))
## `summarise()` has grouped output by 'spp_name_phylo'. You can override using the `.groups` argument.
      # remove underscores species names
      pers_bird$spp_name_phylo = gsub("_", " ", pers_bird$spp_name_phylo)
      
      # create matrix
      pers_bird <- data.frame(pers_bird %>% 
                                spread(personality_trait, n, fill = 0))
      
      # set species name as row name for matrix
      row.names(pers_bird) <- pers_bird$spp_name_phylo 
    
      
      pers_bird <- pers_bird[,2:6]

  # matrix    
    birds_matrix <- data.matrix(pers_bird) 
    
# FINAL TREE
  # making the tree
p_b1 <- ggtree(pruned.birdtree, size = 0.3, layout = 'circular', branch.length = 'none') %<+% bird_data + 
  xlim(-40, NA) + 
  geom_tippoint(aes(color=SSD_index)) + 
  scale_color_gradient2(midpoint = 0, low = "red3", mid = "seashell2", high = "deepskyblue2") + 
  geom_tiplab2(size = 2.2, offset = 4, colour = "black", fontface = "italic") +
  theme(legend.position = 'right')
## Warning: `data_frame()` was deprecated in tibble 1.1.0.
## Please use `tibble()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.
## Warning: `mutate_()` was deprecated in dplyr 0.7.0.
## Please use `mutate()` instead.
## See vignette('programming') for more help
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.
  # adding heatmap of traits
p_b2 <- gheatmap(p_b1, birds_matrix, offset=68, width=2, low = "white", high = "mediumseagreen", color=NULL,
                 colnames=T, colnames_angle = 60, colnames_offset_y = .1, colnames_offset_x = .2) +
  theme(plot.tag = element_text(size = 2, face = "bold"),
        legend.text = element_text(size = 8))

p_b2

# ggsave("./figs/finished figs/birdphylo.tiff", p_b2, width=8, height = 8, units = "in")

13.4 mammals

# subset dataset to include only mammals
 mammal_data <- as.data.frame(figs_data %>%
    filter(taxo_group == "mammal"))

# setting up the basic tree structure
 
  # load tree, set node colours
    mammaltree <- read.tree("./trees/mammal_species.nwk")

  # prune tree to get rid of species we no longer have data for
    pruned.mammaltree <- drop.tip(mammaltree, setdiff(mammaltree$tip.label, mammal_data$spp_name_phylo)) 

  # remove underscores from tip labels
    pruned.mammaltree$tip.label = gsub("_", " ", pruned.mammaltree$tip.label)

  # set rownames for labelling tips
    rownames(mammal_data) <- pruned.mammaltree$tip.label

  # remove underscores from species name from mammal dataset
    mammal_data$spp_name_phylo = gsub("_", " ", mammal_data$spp_name_phylo)

  # set row names
    row.names(mammal_data) <- mammal_data$spp_name_phylo 

# set window
# dev.new(width=8,height=7,noRStudioGD = TRUE)

# make a matrix of effect sizes (n) for each species for each personality trait to add to our plot!
    # subset dataset
    pers_mammal <- as.data.frame(pers_new %>%
    filter(taxo_group == "mammal")) 
  
  # make this a matrix-style dataframe
      pers_mammal <- data.frame(pers_mammal %>%
      group_by(spp_name_phylo, personality_trait) %>%
      summarise(n = n()))
## `summarise()` has grouped output by 'spp_name_phylo'. You can override using the `.groups` argument.
      # remove underscores species names
    pers_mammal$spp_name_phylo = gsub("_", " ", pers_mammal$spp_name_phylo)  
      
    pers_mammal <- data.frame(pers_mammal %>% 
                                spread(personality_trait, n, fill = 0))
    
    row.names(pers_mammal) <- pers_mammal$spp_name_phylo 
    
    pers_mammal <- pers_mammal[,2:6]

  # matrix    
    mammal_matrix <- data.matrix(pers_mammal) 
    
  # making the tree
p_m1 <- ggtree(pruned.mammaltree, size = 0.3, layout = 'circular', branch.length = 'none') %<+% mammal_data + 
  geom_tippoint(aes(color=SSD_index)) + 
  scale_color_gradient2(midpoint = 0, low = "red3", mid = "seashell2", high = "deepskyblue2") + 
  geom_tiplab2(size = 2.5, offset = 2, colour = "black", fontface = "italic") +
  theme(legend.position = 'right')

  # adding heatmap of traits
p_m2 <- gheatmap(p_m1, mammal_matrix, offset=32, width=1.3, low = "white", high = "mediumseagreen", color=NULL,
                 colnames=T, colnames_angle = 60, colnames_offset_y = .1, colnames_offset_x = .2) +
  theme(plot.tag = element_text(size = 9, face = "bold"),
        legend.text = element_text(size = 8))

p_m2

# ggsave("./figs/finished figs/mammalphylo.tiff", p_m2, width=8, height = 7, units = "in")

13.5 reptiles

# subset dataset to include only reptiles
 rep_data <- as.data.frame(figs_data %>%
    filter(taxo_group == "reptilia"))
 
 row.names(rep_data) <- rep_data$spp_name_phylo 
 
# setting up the basic tree structure
 
  # load tree, set node colours
    reptree <- read.tree("./trees/reptile_species.nwk")

  # prune tree to get rid of species we no longer have data for
    pruned.reptree <- drop.tip(reptree, setdiff(reptree$tip.label, rep_data$spp_name_phylo)) 

  # remove underscores from tip labels
    pruned.reptree$tip.label = gsub("_", " ", pruned.reptree$tip.label)

  # set rownames for labelling tips
    rownames(rep_data) <- pruned.reptree$tip.label

  # remove underscores from species name from mammal dataset
    rep_data$spp_name_phylo = gsub("_", " ", rep_data$spp_name_phylo)

# set window size
#dev.new(width=7,height=5,noRStudioGD = TRUE)

 # tree structure 
p3 <- ggtree(pruned.reptree, branch.length='none', size = 0.3, layout='circular') %<+% rep_data + 
  geom_tippoint(aes(color=SSD_index)) + 
  scale_color_gradient2(midpoint = 0, low = "red3", mid = "seashell2", high = "deepskyblue2") + 
  geom_tiplab2(align=T, linetype=NA, size=2.5, offset=4, hjust=0, colour = "black", fontface = "italic") 

# make a matrix of effect sizes (n) for each species for each personality trait to add to our plot!
  # subset dataset
      pers_rep <- as.data.frame(pers_new %>%
      filter(taxo_group == "reptilia")) 
  
  # make this a matrix-style dataframe
      pers_rep <- data.frame(pers_rep %>%
      group_by(spp_name_phylo, personality_trait) %>%
      summarise(n = n()))
## `summarise()` has grouped output by 'spp_name_phylo'. You can override using the `.groups` argument.
  # remove underscores from species name from mammal dataset
      pers_rep$spp_name_phylo = gsub("_", " ", pers_rep$spp_name_phylo)
          
      pers_rep <- data.frame(pers_rep %>% 
                                spread(personality_trait, n, fill = 0))
      
    row.names(pers_rep) <- pers_rep$spp_name_phylo 
    
    pers_rep <- pers_rep[,2:6]

  # matrix    
    rep_matrix <- data.matrix(pers_rep) 
    
  # add the heatmap data to our plot
rep_plot <- gheatmap(p3, rep_matrix, offset = 40, width = 3.5,
            low = "white", high = "mediumseagreen", color=NULL, 
            colnames_position="top", 
            colnames_angle=60, colnames_offset_y = 0, 
             hjust=0, font.size=3) #just not aligning properly 

rep_plot

# ggsave("./figs/finished figs/repphylo.tiff", rep_plot, width=7, height = 5, units = "in")

13.6 fish

# subset dataset to include only fish
 fish_data <- as.data.frame(figs_data %>%
    filter(taxo_group == "fish"))
 
 # window size
 #dev.new(width=8,height=6,noRStudioGD = TRUE)
 
# setting up the basic tree structure
 
  # load tree
fishtree <- read.tree("./trees/fish_species.nwk")

# prune tree to get rid of species we no longer have data for
  pruned.fishtree <- drop.tip(fishtree, setdiff(fishtree$tip.label, fish_data$spp_name_phylo)) 

# remove underscores from tip labels
  pruned.fishtree$tip.label = gsub("_", " ", pruned.fishtree$tip.label)

# set rownames for labelling tips
  rownames(fish_data) <- pruned.fishtree$tip.label

# remove underscores from species name from fish dataset
  fish_data$spp_name_phylo = gsub("_", " ", fish_data$spp_name_phylo)

  row.names(fish_data) <- fish_data$spp_name_phylo 

# make a matrix of effect sizes (n) for each species for each personality trait to add to our plot!
    # subset dataset
  pers_fish <- as.data.frame(pers_new %>%
  filter(taxo_group == "fish")) 
  
  # make this a matrix-style dataframe
      pers_fish <- data.frame(pers_fish %>%
      group_by(spp_name_phylo, personality_trait) %>%
      summarise(n = n()))
## `summarise()` has grouped output by 'spp_name_phylo'. You can override using the `.groups` argument.
      # remove underscores from tip labels
    pers_fish$spp_name_phylo = gsub("_", " ", pers_fish$spp_name_phylo)
      
      pers_fish <- data.frame(pers_fish %>% 
                                spread(personality_trait, n, fill = 0))
    
      row.names(pers_fish) <- pers_fish$spp_name_phylo 
    
    pers_fish <- pers_fish[,2:6]

  # matrix    
    fish_matrix <- data.matrix(pers_fish) 
    
    # FINAL TREE   
  p_f1 <- ggtree(pruned.fishtree, size = 0.3, layout = 'circular', branch.length = 'none') %<+% fish_data + 
  xlim(-30, NA) + 
  geom_tippoint(aes(color=SSD_index)) + 
  scale_color_gradient2(midpoint = 0, low = "red3", mid = "seashell2", high = "deepskyblue2") + 
  geom_tiplab2(size = 2.5, offset = 6, colour = "black", fontface = "italic") +
  theme(legend.position = 'right')
  
  # add the heatmap data to our plot
     fish_plot2 <- gheatmap(p_f1, fish_matrix, offset = 170, width = 5.5,
         low = "white", high = "mediumseagreen", color=NULL, 
         colnames_position="bottom", 
         colnames_angle=60, 
         hjust=0, font.size=3) 

fish_plot2     

# ggsave("./figs/finished figs/fishphylo.tiff", fish_plot2, width=8, height = 6, units = "in")

13.7 inverts

# subset dataset to include only inverts
  invert_data <- as.data.frame(figs_data %>%
    filter(taxo_group == "invertebrate"))
 
 # setting up the basic tree structure
 
  # load tree, set node colours
  inverttree <- read.tree("./trees/invert_species.nwk")

  # prune tree to get rid of species we no longer have data for
  pruned.inverttree <- drop.tip(inverttree, setdiff(inverttree$tip.label, invert_data$spp_name_phylo)) 

  # remove underscores from tip labels
  pruned.inverttree$tip.label = gsub("_", " ", pruned.inverttree$tip.label)

  # remove underscores from dataset and fix row names
  invert_data$spp_name_phylo = gsub("_", " ", invert_data$spp_name_phylo)

  row.names(invert_data) <- invert_data$spp_name_phylo 

# set rownames for labelling tips
  rownames(invert_data) <- pruned.inverttree$tip.label
  
# dev.new(width=8,height=6,noRStudioGD = TRUE) 

 # tree structure (cladogram, circular)
  p5 <- ggtree(pruned.inverttree, branch.length='none', layout='circular') %<+% invert_data + 
  geom_tippoint(aes(color=SSD_index)) + 
  scale_color_gradient2(midpoint = 0, low = "red3", mid = "seashell2", high = "deepskyblue2") + 
  geom_tiplab2(align=T, linetype=NA, size=2.2, offset=2, fontface = "italic") +
  theme(legend.position = "right")

# make a matrix of effect sizes (n) for each species for each personality trait to add to our plot!
    # subset dataset
  pers_invert <- as.data.frame(pers_new %>%
  filter(taxo_group == "invertebrate")) 
  
  # make this a matrix-style dataframe
      pers_invert <- data.frame(pers_invert %>%
      group_by(spp_name_phylo, personality_trait) %>%
      summarise(n = n()))
## `summarise()` has grouped output by 'spp_name_phylo'. You can override using the `.groups` argument.
  # remove underscores from dataset and fix row names
      pers_invert$spp_name_phylo = gsub("_", " ", pers_invert$spp_name_phylo)
       
      pers_invert <- data.frame(pers_invert %>% 
                                spread(personality_trait, n, fill = 0))
      
      row.names(pers_invert) <- pers_invert$spp_name_phylo 
    
      pers_invert <- pers_invert[,2:6]

  # matrix    
    invert_matrix <- data.matrix(pers_invert) 
    
  # add the heatmap data to our plot
     invertplot <- gheatmap(p5, invert_matrix, offset = 40, width = 1.5, 
         low = "white", high = "mediumseagreen", color=NULL, 
         colnames_position="bottom", 
         colnames_angle=45, colnames_offset_y = 0, 
         hjust=0, font.size=2.5)

invertplot     

# save plot
# ggsave("./figs/finished figs/invertphylo.tiff", invertplot, width=8, height = 6, units = "in")

These plots were edited together outside of R with the addition of creative commons animal silhouettes from PhyloPic to create Figures 2-6. Figure 1, the PRISMA diagram, was created using sankeymatic.com